"""
Distributions to use in ARCH models. All distributions must inherit from
:class:`Distribution` and provide the same methods with the same inputs.
"""
from __future__ import annotations
from abc import ABCMeta, abstractmethod
from collections.abc import Sequence
from typing import Callable
import warnings
from numpy import (
abs,
array,
asarray,
empty,
exp,
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
from arch.utility.array import AbstractDocStringInheritor, ensure1d
__all__ = ["Distribution", "Normal", "StudentsT", "SkewStudent", "GeneralizedError"]
[docs]
class Distribution(metaclass=ABCMeta):
"""
Template for subclassing only
"""
def __init__(
self,
random_state: RandomState | None = None,
*,
seed: None | int | RandomState | Generator = None,
) -> None:
self._name = "Distribution"
self.num_params: int = 0
self._parameters: Float64Array | None = 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: None | RandomState | Generator | int = random_state
else:
_seed = seed
if _seed is None:
self._generator: 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: None | Sequence[float] | ArrayLike1D
) -> Float64Array:
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(
"{} does not satisfy the bounds requirement "
"of ({}, {})".format(n, *b)
)
return asarray(params)
@property
def generator(self) -> RandomState | Generator:
"""The NumPy Generator or RandomState attached to the distribution"""
return self._generator
@property
def random_state(self) -> 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: 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: int | float | Sequence[float | int] | ArrayLike1D
) -> Callable[[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: Float64Array) -> 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: Sequence[float] | ArrayLike1D,
resids: ArrayLike,
sigma2: ArrayLike,
individual: bool = False,
) -> float | Float64Array:
"""
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: Float64Array) -> Float64Array:
"""
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: float | Sequence[float] | ArrayLike1D,
parameters: None | Sequence[float] | ArrayLike1D = None,
) -> 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: Sequence[float] | ArrayLike1D,
parameters: None | Sequence[float] | ArrayLike1D = 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: None | Sequence[float] | ArrayLike1D = 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: None | Sequence[float] | ArrayLike1D = 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: RandomState | None = None,
*,
seed: 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: Float64Array) -> list[tuple[float, float]]:
return []
[docs]
def loglikelihood(
self,
parameters: Sequence[float] | ArrayLike1D,
resids: ArrayLike,
sigma2: ArrayLike,
individual: bool = False,
) -> float | Float64Array:
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: Float64Array) -> Float64Array:
return empty(0)
def _simulator(self, size: int | tuple[int, ...]) -> Float64Array:
return self._generator.standard_normal(size)
[docs]
def simulate(
self, parameters: int | float | Sequence[float | int] | ArrayLike1D
) -> Callable[[int | tuple[int, ...]], Float64Array]:
return self._simulator
[docs]
def parameter_names(self) -> list[str]:
return []
[docs]
def cdf(
self,
resids: Sequence[float] | ArrayLike1D,
parameters: None | Sequence[float] | ArrayLike1D = None,
) -> Float64Array:
self._check_constraints(parameters)
return stats.norm.cdf(asarray(resids))
[docs]
def ppf(
self,
pits: float | Sequence[float] | ArrayLike1D,
parameters: None | Sequence[float] | ArrayLike1D = None,
) -> Float64Array:
self._check_constraints(parameters)
scalar = isscalar(pits)
if scalar:
pits = array([pits])
else:
pits = asarray(pits)
ppf = stats.norm.ppf(pits)
if scalar:
return ppf[0]
else:
return ppf
[docs]
def moment(
self, n: int, parameters: None | Sequence[float] | ArrayLike1D = None
) -> float:
if n < 0:
return nan
return stats.norm.moment(n)
[docs]
def partial_moment(
self,
n: int,
z: float = 0.0,
parameters: None | Sequence[float] | ArrayLike1D = 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: RandomState | None = None,
*,
seed: 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: Float64Array) -> list[tuple[float, float]]:
return [(2.05, 500.0)]
[docs]
def loglikelihood(
self,
parameters: Sequence[float] | ArrayLike1D,
resids: ArrayLike,
sigma2: ArrayLike,
individual: bool = False,
) -> float | Float64Array:
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: Float64Array) -> Float64Array:
"""
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: 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: int | float | Sequence[float | int] | ArrayLike1D
) -> Callable[[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 = asarray(parameters, dtype=float)
return self._simulator
[docs]
def parameter_names(self) -> list[str]:
return ["nu"]
[docs]
def cdf(
self,
resids: Sequence[float] | ArrayLike1D,
parameters: None | Sequence[float] | ArrayLike1D = 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: float | Sequence[float] | ArrayLike1D,
parameters: None | Sequence[float] | ArrayLike1D = 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: None | Sequence[float] | ArrayLike1D = 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: None | Sequence[float] | ArrayLike1D = 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: RandomState | None = None,
*,
seed: 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: Float64Array) -> list[tuple[float, float]]:
return [(2.05, 300.0), (-1, 1)]
[docs]
def loglikelihood(
self,
parameters: Sequence[float] | ArrayLike1D,
resids: ArrayLike,
sigma2: ArrayLike,
individual: bool = False,
) -> Float64Array:
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 = asarray(parameters, dtype=float)
eta, lam = parameters
const_c = self.__const_c(parameters)
const_a = self.__const_a(parameters)
const_b = self.__const_b(parameters)
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: Float64Array) -> Float64Array:
"""
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: 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: int | float | Sequence[float | int] | ArrayLike1D
) -> Callable[[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 = asarray(parameters, dtype=float)
return self._simulator
[docs]
def parameter_names(self) -> list[str]:
return ["eta", "lambda"]
def __const_a(self, parameters: Float64Array | Sequence[float]) -> 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: Float64Array | Sequence[float]) -> 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: Float64Array | Sequence[float]) -> 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: Sequence[float] | ArrayLike1D,
parameters: None | Sequence[float] | ArrayLike1D = None,
) -> Float64Array:
parameters = self._check_constraints(parameters)
scalar = isscalar(resids)
if scalar:
resids = array([resids])
resids = asarray(resids, dtype=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
resids = asarray(resids)
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: float | Sequence[float] | ArrayLike1D,
parameters: None | Sequence[float] | ArrayLike1D = None,
) -> 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 / b
if scalar:
return float(icdf[0])
assert isinstance(icdf, ndarray)
return icdf
[docs]
def moment(
self, n: int, parameters: None | Sequence[float] | ArrayLike1D = 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: None | Sequence[float] | ArrayLike1D = 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
lbound = min(z, loc)
lhs = (
(1 - lam)
* (loc ** (n - k))
* (lscale ** k)
* StudentsT._ord_t_partial_moment(k, z=(lbound - loc) / lscale, nu=eta)
)
if z > loc:
rhs = (
(1 + lam)
* (loc ** (n - k))
* (rscale ** k)
* (
StudentsT._ord_t_partial_moment(k, z=(z - loc) / rscale, nu=eta)
- StudentsT._ord_t_partial_moment(k, z=0.0, nu=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: RandomState | None = None,
*,
seed: 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: Float64Array) -> list[tuple[float, float]]:
return [(1.01, 500.0)]
[docs]
def loglikelihood(
self,
parameters: Sequence[float] | ArrayLike1D,
resids: ArrayLike,
sigma2: ArrayLike,
individual: bool = False,
) -> Float64Array:
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: Float64Array) -> Float64Array:
"""
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 array([1.5])
def _simulator(self, size: 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: int | float | Sequence[float | int] | ArrayLike1D
) -> Callable[[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 = asarray(parameters, dtype=float)
return self._simulator
[docs]
def parameter_names(self) -> list[str]:
return ["nu"]
[docs]
def ppf(
self,
pits: float | Sequence[float] | ArrayLike1D,
parameters: None | Sequence[float] | ArrayLike1D = 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: Sequence[float] | ArrayLike1D,
parameters: None | Sequence[float] | ArrayLike1D = None,
) -> Float64Array:
parameters = self._check_constraints(parameters)
nu = parameters[0]
var = stats.gennorm(nu).var()
resids = asarray(resids, dtype=float)
return stats.gennorm(nu, scale=1.0 / sqrt(var)).cdf(resids)
[docs]
def moment(
self, n: int, parameters: None | Sequence[float] | ArrayLike1D = 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: None | Sequence[float] | ArrayLike1D = 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