"""
Asymmetric kernels for R+ and unit interval
References
----------
.. [1] Bouezmarni, Taoufik, and Olivier Scaillet. 2005. “Consistency of
Asymmetric Kernel Density Estimators and Smoothed Histograms with
Application to Income Data.” Econometric Theory 21 (2): 390-412.
.. [2] Chen, Song Xi. 1999. “Beta Kernel Estimators for Density Functions.”
Computational Statistics & Data Analysis 31 (2): 131-45.
https://doi.org/10.1016/S0167-9473(99)00010-9.
.. [3] Chen, Song Xi. 2000. “Probability Density Function Estimation Using
Gamma Kernels.”
Annals of the Institute of Statistical Mathematics 52 (3): 471-80.
https://doi.org/10.1023/A:1004165218295.
.. [4] Jin, Xiaodong, and Janusz Kawczak. 2003. “Birnbaum-Saunders and
Lognormal Kernel Estimators for Modelling Durations in High Frequency
Financial Data.” Annals of Economics and Finance 4: 103-24.
.. [5] Micheaux, Pierre Lafaye de, and Frédéric Ouimet. 2020. “A Study of Seven
Asymmetric Kernels for the Estimation of Cumulative Distribution Functions,”
November. https://arxiv.org/abs/2011.14893v1.
.. [6] Mombeni, Habib Allah, B Masouri, and Mohammad Reza Akhoond. 2019.
“Asymmetric Kernels for Boundary Modification in Distribution Function
Estimation.” REVSTAT, 1-27.
.. [7] Scaillet, O. 2004. “Density Estimation Using Inverse and Reciprocal
Inverse Gaussian Kernels.”
Journal of Nonparametric Statistics 16 (1-2): 217-26.
https://doi.org/10.1080/10485250310001624819.
Created on Mon Mar 8 11:12:24 2021
Author: Josef Perktold
License: BSD-3
"""
import numpy as np
from scipy import special, stats
from statsmodels.tools.validation import float_like
doc_params = """\
Parameters
----------
x : float or array_like
Points at which the kernel is evaluated. ``x`` can be scalar or
1-dim.
sample : array_like
1-d sample from which the kernel estimate is computed.
bw : float
Bandwidth parameter, there is currently no default value for it.
Returns
-------
ndarray
Kernel values evaluated at `x` for each point in `sample`."""
[docs]
def pdf_kernel_asym(x, sample, bw, kernel_type, weights=None, batch_size=10):
"""Density estimate based on asymmetric kernel
Parameters
----------
x : float or array_like
Points for which density is evaluated. ``x`` can be scalar or 1-dim.
sample : array_like
1-d sample from which kernel estimate is computed.
bw : float
Bandwidth parameter, there is currently no default value for it.
kernel_type : str or callable
Kernel name or kernel function.
Currently supported kernel names are "beta", "beta2", "gamma",
"gamma2", "bs", "invgamma", "invgauss", "lognorm", "recipinvgauss" and
"weibull".
weights : None or array_like, optional
If weights is not None, then kernel for sample points are weighted
by it. No weights corresponds to uniform weighting of each component
with 1 / nobs, where nobs is the size of `sample`.
batch_size : float, optional
If x is an 1-dim array, then points can be evaluated in vectorized
form. To limit the amount of memory, a loop works in batches.
Evaluation points and sample points are split so that the
intermediate kernel arrays have at most ``batch_size * 1000``
elements. The value does not need to be an integer.
Default is to have at most 10000 elements in intermediate arrays.
.. versionchanged:: 0.15.1
``batch_size`` bounds the size of the intermediate kernel arrays
on both axes, the evaluation points and the sample. Previously a
sample with more than ``batch_size * 1000`` observations raised an
error. A single-element ``x`` now always returns a scalar.
Returns
-------
pdf : float or ndarray
Estimate of pdf at points x. A float is returned if ``x`` is a scalar
or has a single element, otherwise an ndarray with one entry for each
element of ``x``.
Raises
------
TypeError
If ``batch_size`` is not a number, for example None, a str or a bool.
ValueError
If ``batch_size * 1000`` is smaller than 1 or is nan.
"""
if callable(kernel_type):
kfunc = kernel_type
else:
kfunc = kernel_dict_pdf[kernel_type]
max_elements = _check_batch_size(batch_size)
if np.size(x) * len(sample) < max_elements:
# no batch-loop
if np.size(x) > 1:
x = np.asarray(x)[:, None]
elif np.ndim(x) > 0:
# one-element list, tuple, Series or array: an array is needed
# for the arithmetic inside the kernel. Scalars and 0-dim arrays
# are passed on unchanged.
x = np.asarray(x)
pdfi = kfunc(x, sample, bw)
if weights is None:
pdf = pdfi.mean(-1)
else:
pdf = pdfi @ weights
else:
pdf = _kernel_asym_batch(
kfunc, np.atleast_1d(x), sample, bw, weights, max_elements
)
return pdf
[docs]
def cdf_kernel_asym(x, sample, bw, kernel_type, weights=None, batch_size=10):
"""Estimate of cumulative distribution based on asymmetric kernel
Parameters
----------
x : float or array_like
Points for which cdf is evaluated. ``x`` can be scalar or 1-dim.
sample : array_like
1-d sample from which kernel estimate is computed.
bw : float
Bandwidth parameter, there is currently no default value for it.
kernel_type : str or callable
Kernel name or kernel function.
Currently supported kernel names are "beta", "beta2", "gamma",
"gamma2", "bs", "invgamma", "invgauss", "lognorm", "recipinvgauss" and
"weibull".
weights : None or array_like, optional
If weights is not None, then kernel for sample points are weighted
by it. No weights corresponds to uniform weighting of each component
with 1 / nobs, where nobs is the size of `sample`.
batch_size : float, optional
If x is an 1-dim array, then points can be evaluated in vectorized
form. To limit the amount of memory, a loop works in batches.
Evaluation points and sample points are split so that the
intermediate kernel arrays have at most ``batch_size * 1000``
elements. The value does not need to be an integer.
Default is to have at most 10000 elements in intermediate arrays.
.. versionchanged:: 0.15.1
``batch_size`` bounds the size of the intermediate kernel arrays
on both axes, the evaluation points and the sample. Previously a
sample with more than ``batch_size * 1000`` observations raised an
error. A single-element ``x`` now always returns a scalar.
Returns
-------
cdf : float or ndarray
Estimate of cdf at points x. A float is returned if ``x`` is a scalar
or has a single element, otherwise an ndarray with one entry for each
element of ``x``.
Raises
------
TypeError
If ``batch_size`` is not a number, for example None, a str or a bool.
ValueError
If ``batch_size * 1000`` is smaller than 1 or is nan.
"""
if callable(kernel_type):
kfunc = kernel_type
else:
kfunc = kernel_dict_cdf[kernel_type]
max_elements = _check_batch_size(batch_size)
if np.size(x) * len(sample) < max_elements:
# no batch-loop
if np.size(x) > 1:
x = np.asarray(x)[:, None]
elif np.ndim(x) > 0:
# one-element list, tuple, Series or array: an array is needed
# for the arithmetic inside the kernel. Scalars and 0-dim arrays
# are passed on unchanged.
x = np.asarray(x)
cdfi = kfunc(x, sample, bw)
if weights is None:
cdf = cdfi.mean(-1)
else:
cdf = cdfi @ weights
else:
cdf = _kernel_asym_batch(
kfunc, np.atleast_1d(x), sample, bw, weights, max_elements
)
return cdf
def _check_batch_size(batch_size):
"""Return the maximum number of elements of an intermediate kernel array.
Parameters
----------
batch_size : float
Batch size in units of 1000 elements.
Returns
-------
float
``batch_size * 1000``, which is at least 1.
Raises
------
TypeError
If `batch_size` is not float_like, for example None, a str or a bool.
ValueError
If ``batch_size * 1000`` is smaller than 1 or is nan.
"""
batch_size = float_like(batch_size, "batch_size")
max_elements = batch_size * 1000
# written so that nan fails the check
if not max_elements >= 1:
raise ValueError(
"batch_size * 1000 must be at least 1, the maximum number of "
f"elements in an intermediate array. Got batch_size={batch_size}."
)
return max_elements
def _kernel_asym_batch(kfunc, x, sample, bw, weights, max_elements):
"""Evaluate the weighted kernel sum in blocks of bounded size.
Evaluation points and sample points are split so that each call of
``kfunc`` handles at most ``max_elements`` elements. If the sample alone
exceeds that budget, one evaluation point is combined with a sample block
of at most ``max_elements`` points.
Parameters
----------
kfunc : callable
Kernel function returning an array with one row per evaluation point
and one column per sample point.
x : ndarray
Array of evaluation points, typically 1-dim.
sample : array_like
1-dim sample from which the kernel estimate is computed.
bw : float
Bandwidth parameter.
weights : None or array_like
If None, the kernel values are averaged over the sample. Otherwise
they are summed with these weights, without normalizing them.
max_elements : int or float
Upper bound for the number of elements of the intermediate kernel
arrays, that is ``batch_size * 1000``.
Returns
-------
float or ndarray
Estimate at the points in `x`. A scalar is returned if `x` has a
single element and at most one dimension.
"""
nobs = len(sample)
# Evaluation points that can be combined with the full sample. A
# non-integer ``batch_size`` was accepted before, so cast for ``range``.
max_x = max(1, int(max_elements // nobs))
# Sample points that can be combined with ``max_x`` evaluation points.
max_obs = max(1, int(max_elements // max_x))
# The sample only has to be materialized when it is split. Otherwise the
# kernel receives the caller's object (Series, list, ...) unchanged.
sample_array = np.asarray(sample) if nobs > max_obs else None
if weights is not None:
weights = np.asarray(weights)
if weights.ndim == 0 or weights.shape[0] != nobs:
raise ValueError("weights must have one entry per sample point")
estimates = []
for start in range(0, len(x), max_x):
x_block = x[start : start + max_x]
chunk_sums = []
for obs_start in range(0, nobs, max_obs):
obs_stop = obs_start + max_obs
if sample_array is None:
sample_chunk = sample
else:
sample_chunk = sample_array[obs_start:obs_stop]
kernel = kfunc(x_block[:, None], sample_chunk, bw)
if weights is None:
chunk_sums.append(kernel.sum(-1))
else:
chunk_sums.append(kernel @ weights[obs_start:obs_stop])
block = np.sum(chunk_sums, axis=0)
if weights is None:
block = block / nobs
estimates.append(block)
estimate = np.concatenate(estimates)
if np.size(x) == 1 and x.ndim <= 1:
estimate = estimate[0]
return estimate
[docs]
def kernel_pdf_beta(x, sample, bw):
# Beta kernel for density, pdf, estimation
return stats.beta.pdf(sample, x / bw + 1, (1 - x) / bw + 1)
kernel_pdf_beta.__doc__ = f"""\
Beta kernel for density, pdf, estimation
{doc_params}
References
----------
.. [1] Bouezmarni, Taoufik, and Olivier Scaillet. 2005. “Consistency of
Asymmetric Kernel Density Estimators and Smoothed Histograms with
Application to Income Data.” Econometric Theory 21 (2): 390-412.
.. [2] Chen, Song Xi. 1999. “Beta Kernel Estimators for Density Functions.”
Computational Statistics & Data Analysis 31 (2): 131-45.
https://doi.org/10.1016/S0167-9473(99)00010-9.
"""
[docs]
def kernel_cdf_beta(x, sample, bw):
# Beta kernel for cumulative distribution, cdf, estimation
return stats.beta.sf(sample, x / bw + 1, (1 - x) / bw + 1)
kernel_cdf_beta.__doc__ = f"""\
Beta kernel for cumulative distribution, cdf, estimation
{doc_params}
References
----------
.. [1] Bouezmarni, Taoufik, and Olivier Scaillet. 2005. “Consistency of
Asymmetric Kernel Density Estimators and Smoothed Histograms with
Application to Income Data.” Econometric Theory 21 (2): 390-412.
.. [2] Chen, Song Xi. 1999. “Beta Kernel Estimators for Density Functions.”
Computational Statistics & Data Analysis 31 (2): 131-45.
https://doi.org/10.1016/S0167-9473(99)00010-9.
"""
[docs]
def kernel_pdf_beta2(x, sample, bw):
# Beta kernel for density, pdf, estimation with boundary corrections
# a = 2 * bw**2 + 2.5 -
# np.sqrt(4 * bw**4 + 6 * bw**2 + 2.25 - x**2 - x / bw)
# terms a1 and a2 are independent of x
a1 = 2 * bw**2 + 2.5
a2 = 4 * bw**4 + 6 * bw**2 + 2.25
if np.size(x) == 1:
# without vectorizing:
if x < 2 * bw:
a = a1 - np.sqrt(a2 - x**2 - x / bw)
pdf = stats.beta.pdf(sample, a, (1 - x) / bw)
elif x > (1 - 2 * bw):
x_ = 1 - x
a = a1 - np.sqrt(a2 - x_**2 - x_ / bw)
pdf = stats.beta.pdf(sample, x / bw, a)
else:
pdf = stats.beta.pdf(sample, x / bw, (1 - x) / bw)
else:
alpha = x / bw
beta = (1 - x) / bw
mask_low = x < 2 * bw
x_ = x[mask_low]
alpha[mask_low] = a1 - np.sqrt(a2 - x_**2 - x_ / bw)
mask_upp = x > (1 - 2 * bw)
x_ = 1 - x[mask_upp]
beta[mask_upp] = a1 - np.sqrt(a2 - x_**2 - x_ / bw)
pdf = stats.beta.pdf(sample, alpha, beta)
return pdf
kernel_pdf_beta2.__doc__ = f"""\
Beta kernel for density, pdf, estimation with boundary corrections
{doc_params}
References
----------
.. [1] Bouezmarni, Taoufik, and Olivier Scaillet. 2005. “Consistency of
Asymmetric Kernel Density Estimators and Smoothed Histograms with
Application to Income Data.” Econometric Theory 21 (2): 390-412.
.. [2] Chen, Song Xi. 1999. “Beta Kernel Estimators for Density Functions.”
Computational Statistics & Data Analysis 31 (2): 131-45.
https://doi.org/10.1016/S0167-9473(99)00010-9.
"""
[docs]
def kernel_cdf_beta2(x, sample, bw):
# Beta kernel for cdf estimation with boundary correction
# a = 2 * bw**2 + 2.5 -
# np.sqrt(4 * bw**4 + 6 * bw**2 + 2.25 - x**2 - x / bw)
# terms a1 and a2 are independent of x
a1 = 2 * bw**2 + 2.5
a2 = 4 * bw**4 + 6 * bw**2 + 2.25
if np.size(x) == 1:
# without vectorizing:
if x < 2 * bw:
a = a1 - np.sqrt(a2 - x**2 - x / bw)
pdf = stats.beta.sf(sample, a, (1 - x) / bw)
elif x > (1 - 2 * bw):
x_ = 1 - x
a = a1 - np.sqrt(a2 - x_**2 - x_ / bw)
pdf = stats.beta.sf(sample, x / bw, a)
else:
pdf = stats.beta.sf(sample, x / bw, (1 - x) / bw)
else:
alpha = x / bw
beta = (1 - x) / bw
mask_low = x < 2 * bw
x_ = x[mask_low]
alpha[mask_low] = a1 - np.sqrt(a2 - x_**2 - x_ / bw)
mask_upp = x > (1 - 2 * bw)
x_ = 1 - x[mask_upp]
beta[mask_upp] = a1 - np.sqrt(a2 - x_**2 - x_ / bw)
pdf = stats.beta.sf(sample, alpha, beta)
return pdf
kernel_cdf_beta2.__doc__ = f"""\
Beta kernel for cdf estimation with boundary correction
{doc_params}
References
----------
.. [1] Bouezmarni, Taoufik, and Olivier Scaillet. 2005. “Consistency of
Asymmetric Kernel Density Estimators and Smoothed Histograms with
Application to Income Data.” Econometric Theory 21 (2): 390-412.
.. [2] Chen, Song Xi. 1999. “Beta Kernel Estimators for Density Functions.”
Computational Statistics & Data Analysis 31 (2): 131-45.
https://doi.org/10.1016/S0167-9473(99)00010-9.
"""
[docs]
def kernel_pdf_gamma(x, sample, bw):
# Gamma kernel for density, pdf, estimation
pdfi = stats.gamma.pdf(sample, x / bw + 1, scale=bw)
return pdfi
kernel_pdf_gamma.__doc__ = f"""\
Gamma kernel for density, pdf, estimation
{doc_params}
References
----------
.. [1] Bouezmarni, Taoufik, and Olivier Scaillet. 2005. “Consistency of
Asymmetric Kernel Density Estimators and Smoothed Histograms with
Application to Income Data.” Econometric Theory 21 (2): 390-412.
.. [2] Chen, Song Xi. 2000. “Probability Density Function Estimation Using
Gamma Kernels.”
Annals of the Institute of Statistical Mathematics 52 (3): 471-80.
https://doi.org/10.1023/A:1004165218295.
"""
[docs]
def kernel_cdf_gamma(x, sample, bw):
# Gamma kernel for density, pdf, estimation
# kernel cdf uses the survival function, but I don't know why.
cdfi = stats.gamma.sf(sample, x / bw + 1, scale=bw)
return cdfi
kernel_cdf_gamma.__doc__ = f"""\
Gamma kernel for cumulative distribution, cdf, estimation
{doc_params}
References
----------
.. [1] Bouezmarni, Taoufik, and Olivier Scaillet. 2005. “Consistency of
Asymmetric Kernel Density Estimators and Smoothed Histograms with
Application to Income Data.” Econometric Theory 21 (2): 390-412.
.. [2] Chen, Song Xi. 2000. “Probability Density Function Estimation Using
Gamma Kernels.”
Annals of the Institute of Statistical Mathematics 52 (3): 471-80.
https://doi.org/10.1023/A:1004165218295.
"""
def _kernel_pdf_gamma(x, sample, bw):
"""Gamma kernel for pdf, without boundary corrected part
drops `+ 1` in shape parameter
It should be possible to use this if probability in
neighborhood of zero boundary is small.
Parameters
----------
x : float or array_like
Points for which density is evaluated. ``x`` can be scalar or 1-dim.
sample : array_like
1-d sample from which kernel estimate is computed.
bw : float
Bandwidth parameter, there is currently no default value for it.
Returns
-------
pdf : float or ndarray
Estimate of pdf at points x.
"""
return stats.gamma.pdf(sample, x / bw, scale=bw)
def _kernel_cdf_gamma(x, sample, bw):
"""Gamma kernel for cdf, without boundary corrected part
drops `+ 1` in shape parameter
It should be possible to use this if probability in
neighborhood of zero boundary is small.
Parameters
----------
x : float or array_like
Points for which density is evaluated. ``x`` can be scalar or 1-dim.
sample : array_like
1-d sample from which kernel estimate is computed.
bw : float
Bandwidth parameter, there is currently no default value for it.
Returns
-------
cdf : float or ndarray
Estimate of cdf at points x.
"""
return stats.gamma.sf(sample, x / bw, scale=bw)
[docs]
def kernel_pdf_gamma2(x, sample, bw):
# Gamma kernel for density, pdf, estimation with boundary correction
if np.size(x) == 1:
# without vectorizing, easier to read
if x < 2 * bw:
a = (x / bw) ** 2 + 1
else:
a = x / bw
else:
a = x / bw
mask = x < 2 * bw
a[mask] = a[mask] ** 2 + 1
pdf = stats.gamma.pdf(sample, a, scale=bw)
return pdf
kernel_pdf_gamma2.__doc__ = f"""\
Gamma kernel for density, pdf, estimation with boundary correction
{doc_params}
References
----------
.. [1] Bouezmarni, Taoufik, and Olivier Scaillet. 2005. “Consistency of
Asymmetric Kernel Density Estimators and Smoothed Histograms with
Application to Income Data.” Econometric Theory 21 (2): 390-412.
.. [2] Chen, Song Xi. 2000. “Probability Density Function Estimation Using
Gamma Kernels.”
Annals of the Institute of Statistical Mathematics 52 (3): 471-80.
https://doi.org/10.1023/A:1004165218295.
"""
[docs]
def kernel_cdf_gamma2(x, sample, bw):
# Gamma kernel for cdf estimation with boundary correction
if np.size(x) == 1:
# without vectorizing
if x < 2 * bw:
a = (x / bw) ** 2 + 1
else:
a = x / bw
else:
a = x / bw
mask = x < 2 * bw
a[mask] = a[mask] ** 2 + 1
pdf = stats.gamma.sf(sample, a, scale=bw)
return pdf
kernel_cdf_gamma2.__doc__ = f"""\
Gamma kernel for cdf estimation with boundary correction
{doc_params}
References
----------
.. [1] Bouezmarni, Taoufik, and Olivier Scaillet. 2005. “Consistency of
Asymmetric Kernel Density Estimators and Smoothed Histograms with
Application to Income Data.” Econometric Theory 21 (2): 390-412.
.. [2] Chen, Song Xi. 2000. “Probability Density Function Estimation Using
Gamma Kernels.”
Annals of the Institute of Statistical Mathematics 52 (3): 471-80.
https://doi.org/10.1023/A:1004165218295.
"""
[docs]
def kernel_pdf_invgamma(x, sample, bw):
# Inverse gamma kernel for density, pdf, estimation
return stats.invgamma.pdf(sample, 1 / bw + 1, scale=x / bw)
kernel_pdf_invgamma.__doc__ = f"""\
Inverse gamma kernel for density, pdf, estimation
Based on cdf kernel by Micheaux and Ouimet (2020)
{doc_params}
References
----------
.. [1] Micheaux, Pierre Lafaye de, and Frédéric Ouimet. 2020. “A Study of
Seven Asymmetric Kernels for the Estimation of Cumulative Distribution
Functions,” November. https://arxiv.org/abs/2011.14893v1.
"""
[docs]
def kernel_cdf_invgamma(x, sample, bw):
# Inverse gamma kernel for cumulative distribution, cdf, estimation
return stats.invgamma.sf(sample, 1 / bw + 1, scale=x / bw)
kernel_cdf_invgamma.__doc__ = f"""\
Inverse gamma kernel for cumulative distribution, cdf, estimation
{doc_params}
References
----------
.. [1] Micheaux, Pierre Lafaye de, and Frédéric Ouimet. 2020. “A Study of
Seven Asymmetric Kernels for the Estimation of Cumulative Distribution
Functions,” November. https://arxiv.org/abs/2011.14893v1.
"""
[docs]
def kernel_pdf_invgauss(x, sample, bw):
# Inverse gaussian kernel for density, pdf, estimation
m = x
lam = 1 / bw
return stats.invgauss.pdf(sample, m / lam, scale=lam)
kernel_pdf_invgauss.__doc__ = f"""\
Inverse gaussian kernel for density, pdf, estimation
{doc_params}
References
----------
.. [1] Scaillet, O. 2004. “Density Estimation Using Inverse and Reciprocal
Inverse Gaussian Kernels.”
Journal of Nonparametric Statistics 16 (1-2): 217-26.
https://doi.org/10.1080/10485250310001624819.
"""
def kernel_pdf_invgauss_(x, sample, bw):
"""Inverse gaussian kernel density, explicit formula
Scaillet 2004
Parameters
----------
x : float or array_like
Points for which density is evaluated. ``x`` can be scalar or 1-dim.
sample : array_like
1-d sample from which kernel estimate is computed.
bw : float
Bandwidth parameter, there is currently no default value for it.
Returns
-------
pdf : float or ndarray
Estimate of pdf at points x, averaged over `sample`.
"""
pdf = (
1
/ np.sqrt(2 * np.pi * bw * sample**3)
* np.exp(-1 / (2 * bw * x) * (sample / x - 2 + x / sample))
)
return pdf.mean(-1)
[docs]
def kernel_cdf_invgauss(x, sample, bw):
# Inverse gaussian kernel for cumulative distribution, cdf, estimation
m = x
lam = 1 / bw
return stats.invgauss.sf(sample, m / lam, scale=lam)
kernel_cdf_invgauss.__doc__ = f"""\
Inverse gaussian kernel for cumulative distribution, cdf, estimation
{doc_params}
References
----------
.. [1] Scaillet, O. 2004. “Density Estimation Using Inverse and Reciprocal
Inverse Gaussian Kernels.”
Journal of Nonparametric Statistics 16 (1-2): 217-26.
https://doi.org/10.1080/10485250310001624819.
"""
[docs]
def kernel_pdf_recipinvgauss(x, sample, bw):
# Reciprocal inverse gaussian kernel for density, pdf, estimation
# need shape-scale parameterization for scipy
# references use m, lambda parameterization
m = 1 / (x - bw)
lam = 1 / bw
return stats.recipinvgauss.pdf(sample, m / lam, scale=1 / lam)
kernel_pdf_recipinvgauss.__doc__ = f"""\
Reciprocal inverse gaussian kernel for density, pdf, estimation
{doc_params}
References
----------
.. [1] Scaillet, O. 2004. “Density Estimation Using Inverse and Reciprocal
Inverse Gaussian Kernels.”
Journal of Nonparametric Statistics 16 (1-2): 217-26.
https://doi.org/10.1080/10485250310001624819.
"""
def kernel_pdf_recipinvgauss_(x, sample, bw):
"""Reciprocal inverse gaussian kernel density, explicit formula
Scaillet 2004
Parameters
----------
x : float or array_like
Points for which density is evaluated. ``x`` can be scalar or 1-dim.
sample : array_like
1-d sample from which kernel estimate is computed.
bw : float
Bandwidth parameter, there is currently no default value for it.
Returns
-------
pdf : float or ndarray
Kernel component of the pdf estimate at points x, one value for
each point in `sample`.
"""
pdf = (
1
/ np.sqrt(2 * np.pi * bw * sample)
* np.exp(-(x - bw) / (2 * bw) * sample / (x - bw) - 2 + (x - bw) / sample)
)
return pdf
[docs]
def kernel_cdf_recipinvgauss(x, sample, bw):
# Reciprocal inverse gaussian kernel for cdf estimation
# need shape-scale parameterization for scipy
# references use m, lambda parameterization
m = 1 / (x - bw)
lam = 1 / bw
return stats.recipinvgauss.sf(sample, m / lam, scale=1 / lam)
kernel_cdf_recipinvgauss.__doc__ = f"""\
Reciprocal inverse gaussian kernel for cdf estimation
{doc_params}
References
----------
.. [1] Scaillet, O. 2004. “Density Estimation Using Inverse and Reciprocal
Inverse Gaussian Kernels.”
Journal of Nonparametric Statistics 16 (1-2): 217-26.
https://doi.org/10.1080/10485250310001624819.
"""
[docs]
def kernel_pdf_bs(x, sample, bw):
# Birnbaum Saunders (normal) kernel for density, pdf, estimation
return stats.fatiguelife.pdf(sample, bw, scale=x)
kernel_pdf_bs.__doc__ = f"""\
Birnbaum Saunders (normal) kernel for density, pdf, estimation
{doc_params}
References
----------
.. [1] Jin, Xiaodong, and Janusz Kawczak. 2003. “Birnbaum-Saunders and
Lognormal Kernel Estimators for Modelling Durations in High Frequency
Financial Data.” Annals of Economics and Finance 4: 103-24.
"""
[docs]
def kernel_cdf_bs(x, sample, bw):
# Birnbaum Saunders (normal) kernel for cdf estimation
return stats.fatiguelife.sf(sample, bw, scale=x)
kernel_cdf_bs.__doc__ = f"""\
Birnbaum Saunders (normal) kernel for cdf estimation
{doc_params}
References
----------
.. [1] Jin, Xiaodong, and Janusz Kawczak. 2003. “Birnbaum-Saunders and
Lognormal Kernel Estimators for Modelling Durations in High Frequency
Financial Data.” Annals of Economics and Finance 4: 103-24.
.. [2] Mombeni, Habib Allah, B Masouri, and Mohammad Reza Akhoond. 2019.
“Asymmetric Kernels for Boundary Modification in Distribution Function
Estimation.” REVSTAT, 1-27.
"""
[docs]
def kernel_pdf_lognorm(x, sample, bw):
# Log-normal kernel for density, pdf, estimation
# need shape-scale parameterization for scipy
# not sure why JK picked this normalization, makes required bw small
# maybe we should skip this transformation and just use bw
# Funke and Kawka 2015 (table 1) use bw (or bw**2) corresponding to
# variance of normal pdf
# bw = np.exp(bw_**2 / 4) - 1 # this is inverse transformation
bw_ = np.sqrt(4 * np.log(1 + bw))
return stats.lognorm.pdf(sample, bw_, scale=x)
kernel_pdf_lognorm.__doc__ = f"""\
Log-normal kernel for density, pdf, estimation
{doc_params}
Notes
-----
Warning: parameterization of bandwidth will likely be changed
References
----------
.. [1] Jin, Xiaodong, and Janusz Kawczak. 2003. “Birnbaum-Saunders and
Lognormal Kernel Estimators for Modelling Durations in High Frequency
Financial Data.” Annals of Economics and Finance 4: 103-24.
"""
[docs]
def kernel_cdf_lognorm(x, sample, bw):
# Log-normal kernel for cumulative distribution, cdf, estimation
# need shape-scale parameterization for scipy
# not sure why JK picked this normalization, makes required bw small
# maybe we should skip this transformation and just use bw
# Funke and Kawka 2015 (table 1) use bw (or bw**2) corresponding to
# variance of normal pdf
# bw = np.exp(bw_**2 / 4) - 1 # this is inverse transformation
bw_ = np.sqrt(4 * np.log(1 + bw))
return stats.lognorm.sf(sample, bw_, scale=x)
kernel_cdf_lognorm.__doc__ = f"""\
Log-normal kernel for cumulative distribution, cdf, estimation
{doc_params}
Notes
-----
Warning: parameterization of bandwidth will likely be changed
References
----------
.. [1] Jin, Xiaodong, and Janusz Kawczak. 2003. “Birnbaum-Saunders and
Lognormal Kernel Estimators for Modelling Durations in High Frequency
Financial Data.” Annals of Economics and Finance 4: 103-24.
"""
def kernel_pdf_lognorm_(x, sample, bw):
"""Log-normal kernel for density, pdf, estimation, explicit formula
Jin, Kawczak 2003
Parameters
----------
x : float or array_like
Points for which density is evaluated. ``x`` can be scalar or 1-dim.
sample : array_like
1-d sample from which kernel estimate is computed.
bw : float
Bandwidth parameter, there is currently no default value for it.
Returns
-------
pdf : float or ndarray
Estimate of pdf at points x, averaged over `sample`.
"""
term = 8 * np.log(1 + bw) # this is 2 * variance in normal pdf
pdf = (
1
/ np.sqrt(term * np.pi)
/ sample
* np.exp(-((np.log(x) - np.log(sample)) ** 2) / term)
)
return pdf.mean(-1)
[docs]
def kernel_pdf_weibull(x, sample, bw):
# Weibull kernel for density, pdf, estimation
# need shape-scale parameterization for scipy
# references use m, lambda parameterization
return stats.weibull_min.pdf(sample, 1 / bw, scale=x / special.gamma(1 + bw))
kernel_pdf_weibull.__doc__ = f"""\
Weibull kernel for density, pdf, estimation
Based on cdf kernel by Mombeni et al. (2019)
{doc_params}
References
----------
.. [1] Mombeni, Habib Allah, B Masouri, and Mohammad Reza Akhoond. 2019.
“Asymmetric Kernels for Boundary Modification in Distribution Function
Estimation.” REVSTAT, 1-27.
"""
[docs]
def kernel_cdf_weibull(x, sample, bw):
# Weibull kernel for cumulative distribution, cdf, estimation
# need shape-scale parameterization for scipy
# references use m, lambda parameterization
return stats.weibull_min.sf(sample, 1 / bw, scale=x / special.gamma(1 + bw))
kernel_cdf_weibull.__doc__ = f"""\
Weibull kernel for cumulative distribution, cdf, estimation
{doc_params}
References
----------
.. [1] Mombeni, Habib Allah, B Masouri, and Mohammad Reza Akhoond. 2019.
“Asymmetric Kernels for Boundary Modification in Distribution Function
Estimation.” REVSTAT, 1-27.
"""
# produced wth
# print("\n".join(['"%s": %s,' % (i.split("_")[-1], i) for i in dir(kern)
# if "kernel" in i and not i.endswith("_")]))
kernel_dict_cdf = {
"beta": kernel_cdf_beta,
"beta2": kernel_cdf_beta2,
"bs": kernel_cdf_bs,
"gamma": kernel_cdf_gamma,
"gamma2": kernel_cdf_gamma2,
"invgamma": kernel_cdf_invgamma,
"invgauss": kernel_cdf_invgauss,
"lognorm": kernel_cdf_lognorm,
"recipinvgauss": kernel_cdf_recipinvgauss,
"weibull": kernel_cdf_weibull,
}
kernel_dict_pdf = {
"beta": kernel_pdf_beta,
"beta2": kernel_pdf_beta2,
"bs": kernel_pdf_bs,
"gamma": kernel_pdf_gamma,
"gamma2": kernel_pdf_gamma2,
"invgamma": kernel_pdf_invgamma,
"invgauss": kernel_pdf_invgauss,
"lognorm": kernel_pdf_lognorm,
"recipinvgauss": kernel_pdf_recipinvgauss,
"weibull": kernel_pdf_weibull,
}
