cdxcore.rolling#

Robust rolling statistics.

Recursive exponentially weighted estimators of location and scale, in particular robust versions which use the median and MAD for initialization, flag outliers, and do not let them move the location. They are available for regular time steps (cdxcore.rolling.robust_rolling_ew()) and irregular time steps (cdxcore.rolling.robust_rolling_dt_ew()). cdxcore.rolling.robust_fixed_window() is the non-recursive analogue.

import numpy as np
from cdxcore.rolling import robust_rolling_ew, robust_rolling_dt_ew

x = np.random.normal( size=1000 )
x[500] = 10.

loc, vol, outlier = robust_rolling_ew( x, window=20 )          # outlier[500] is True

dt = np.random.uniform( 0.001, 0.01, size=1000 )
dx = 0.1 * dt + 0.2 * np.sqrt( dt ) * np.random.normal( size=1000 )
mu, sigma, outlier = robust_rolling_dt_ew( dx, dt, twindow=0.25, scale_by_dt=True, normalize_by_dt=True )

If numba is installed the recursions are compiled, otherwise they run as plain python loops.

Module Attributes

MAD_NORMAL_SCALE

Gaussian consistency factor for a median absolute deviation.

Functions

`robust_fixed_window`(x[, dt, cutoff, ...])	Robust location and scale of a fixed set of equally weighted samples.
`robust_rolling_dt_ew`(x, dt[, twindow, init, ...])	Computes a robust recursive exponentially weighted mean and volatility for irregular time steps.
`robust_rolling_ew`(x, window[, init, cutoff, ...])	Computes a robust recursive exponentially weighted mean and volatility.
`rolling_ew_std`(x, window[, init, cutoff])	Computes a standard recursive exponentially weighted mean and volatility, initialized over `init` steps.

cdxcore.rolling.MAD_NORMAL_SCALE = 1.482602218505602#: Gaussian consistency factor for a median absolute deviation.

cdxcore.rolling.robust_fixed_window(x, dt=None, cutoff=2.5, scale_by_dt=False, normalize_by_dt=False, min_scale=None)[source]#

Robust location and scale of a fixed set of equally weighted samples.

This is the non-recursive analogue of robust_rolling_dt_ew(): the location is the median, and the scale is the MAD around it, multiplied by 1.4826. There is no warm-up and no outlier mask; cutoff is only retained to match the interface of the recursive functions.

If dt is provided and scale_by_dt is True then the model

x[i] = mu * dt[i] + sigma * sqrt(dt[i]) * eps[i]

is assumed: mu is the median of x/dt, and sigma is the scaled MAD of the residuals (x[i] - mu*dt[i]) / sqrt(dt[i]). If normalize_by_dt is True the function returns mu and sigma. Otherwise it returns mu*dt_ref and sigma*sqrt(dt_ref) where dt_ref is the median of dt, i.e. quantities in the units of x.

Parameters:

xnp.ndarray: One-dimensional samples.
dtnp.ndarray, optional: Time steps, strictly positive. Ignored unless scale_by_dt is True.
cutofffloat, optional: Not used. Default 2.5.
scale_by_dtbool, optional: Use the increment model described above. Requires dt. Default False.
normalize_by_dtbool, optional: Return drift rate and diffusion scale. Requires scale_by_dt. Default False.
min_scalefloat, optional: Floor for the scale. Defaults to machine precision times max(1, max(|x|)).

Returns:

locationfloat: Robust location.
scalefloat: Robust scale, i.e. a standard deviation, not a variance.

cdxcore.rolling.robust_rolling_dt_ew(x, dt, twindow=0.25, init=10, cutoff=2.5, scale_by_dt=False, normalize_by_dt=False, min_scale=None, on_zero_init='warn')[source]#

Computes a robust recursive exponentially weighted mean and volatility for irregular time steps.

The per-step weight is w[t] = 1 - exp(-dt[t]/twindow), which is approximately dt[t]/twindow. With c = 1/E[min(|Z|,cutoff)] for a standard normal Z the update rule is as follows.

If scale_by_dt is False, with the innovation d[t] = x[t] - m[t-1]:

m[t] = m[t-1] + w[t] * d[t]
v[t] = (1-w[t]) * v[t-1] + w[t] * c * abs( clip( d[t], -cutoff*v[t-1], cutoff*v[t-1] ) )

If x is itself a return-type quantity such as dS for a stock, use scale_by_dt=True. The drift m is then a rate, and the innovation is normalized by sqrt(dt):

d[t] = ( x[t] - m[t-1]*dt[t] ) / sqrt(dt[t])
m[t] = m[t-1] + w[t] / dt[t] * ( x[t] - m[t-1]*dt[t] )
v[t] = (1-w[t]) * v[t-1] + w[t] * c * abs( clip( d[t], -cutoff*v[t-1], cutoff*v[t-1] ) )

If all time steps are equal to dt and twindow = window*dt, this is equivalent to robust_rolling_ew() except that the estimated quantity is the mean of dx/dt, and the volatility is that of (dx - m*dt)/sqrt(dt).

An outlier is identified if |d[t]| exceeds cutoff * v[t-1]. In that case m[t] = m[t-1]. The volatility is updated using the capped innovation and is floored at min_scale.

The state at index init-1 is initialized with the median and MAD (scaled by 1.4826) of x[:init], or of x[:init]/dt[:init] and of the normalized innovations if scale_by_dt is True. These are weighted by the weights the observations have in the terminal state of the recursion. Earlier elements are NaN.

Parameters:

xnp.ndarray

One-dimensional time series.

dtnp.ndarray

Time steps, strictly positive, same shape as x.

twindowfloat, optional

Time constant of the weights: a new observation gets the weight 1 - exp(-dt/twindow). Default 0.25.

initint, optional

Initial period. Default 10.

cutofffloat, optional

Normalized innovations exceeding this level are considered outliers. Default 2.5.

scale_by_dtbool, optional

Scale returns by dt and volatilities by sqrt(dt) during estimation, see above. Default False.

normalize_by_dtbool, optional

Requires scale_by_dt; otherwise a ValueError is raised. If True return the drift rate m and volatility v. If False return m*dt and v*sqrt(dt), i.e. quantities in the units of x. Default False.

min_scalefloat, optional

Floor for the volatility, in the units of the normalized innovation. Defaults to machine precision times the largest absolute normalized value of the first init steps (at least 1).

on_zero_initstr, optional

What to do if the median absolute deviation (MAD) of the initial period is zero. This happens if more than half of the (weighted) observations equal the median, e.g. for zero-inflated non-negative series such as variances or volumes. The initial volatility is then only min_scale, any later observation which differs from the location is flagged as an outlier and cannot move it, and the location can stay stuck for a long time. It does not apply if init is 1. Options:

"warn": issue a RuntimeWarning and continue. This is the default.
"ignore": continue silently.
"raise": raise an OverflowError.
"standard": initialize with the (weighted) mean and standard deviation of the initial period instead of the median and MAD. If the initial period is constant the standard deviation is zero as well, and the result is as for "ignore".

Returns:

meannp.ndarray: Robust mean.
volnp.ndarray: Robust volatility, i.e. a standard deviation, not a variance.
outliernp.ndarray: Boolean array which is True where an outlier was detected relative to the previous state. outlier[init-1] is always False.

cdxcore.rolling.robust_rolling_ew(x, window, init=10, cutoff=2.5, min_scale=None, on_zero_init='warn')[source]#

Computes a robust recursive exponentially weighted mean and volatility.

With w = 1/window, the innovation d[t] = x[t] - m[t-1] and c = 1/E[min(|Z|,cutoff)] for a standard normal Z, the update rule is:

m[t] = m[t-1] + w * d[t]
v[t] = (1-w) * v[t-1] + w * c * abs( clip( d[t], -cutoff*v[t-1], cutoff*v[t-1] ) )

An outlier is identified if |d[t]| exceeds cutoff * v[t-1]. In that case the mean is not updated: m[t] = m[t-1]. The volatility is updated using the capped innovation as above and is floored at min_scale.

The state at index init-1 is initialized with the median and MAD (scaled by 1.4826) of x[:init], weighted by the weights the observations have in the terminal state of the recursion. Earlier elements are NaN.

Parameters:

xnp.ndarray

One-dimensional time series.

windowfloat

The parametrization w = 1/window means that any new observation gets the same weight as it would get in a rolling estimator with size window. Must be greater than 1.

initint, optional

Initial period. Default 10.

cutofffloat, optional

Normalized innovations exceeding this level are considered outliers. Default 2.5.

min_scalefloat, optional

Floor for the volatility. Defaults to machine precision times max(1, max(|x[:init]|)).

on_zero_initstr, optional

What to do if the median absolute deviation (MAD) of the initial period is zero. This happens if more than half of the (weighted) observations equal the median, e.g. for zero-inflated non-negative series such as variances or volumes. The initial volatility is then only min_scale, any later observation which differs from the location is flagged as an outlier and cannot move it, and the location can stay stuck for a long time. It does not apply if init is 1. Options:

"warn": issue a RuntimeWarning and continue. This is the default.
"ignore": continue silently.
"raise": raise an OverflowError.
"standard": initialize with the (weighted) mean and standard deviation of the initial period instead of the median and MAD. If the initial period is constant the standard deviation is zero as well, and the result is as for "ignore".

Returns:

meannp.ndarray: Robust mean.
volnp.ndarray: Robust volatility, i.e. a standard deviation, not a variance.
outliernp.ndarray: Boolean array which is True where an outlier was detected.

cdxcore.rolling.rolling_ew_std(x, window, init=10, cutoff=2.5)[source]#

Computes a standard recursive exponentially weighted mean and volatility, initialized over init steps.

With w = 1/window and the innovation d[t] = x[t] - m[t-1] the update rule is:

m[t] = m[t-1] + w * d[t]
v[t] = (1-w) * ( v[t-1] + w * d[t]**2 )

Here v is the variance; the function returns sqrt(v).

An outlier is identified if |d[t]| exceeds cutoff * sqrt(v[t-1]). In that case the mean is not updated, and the variance is updated with the capped innovation:

m[t] = m[t-1]
v[t] = (1-w) * ( v[t-1] + w * clip( d[t], -cutoff*sqrt(v[t-1]), cutoff*sqrt(v[t-1]) )**2 )

The state is initialized at index init-1 with the mean and variance of x[:init]; earlier elements are NaN.

Parameters:

xnp.ndarray: One-dimensional time series.
windowfloat: The parametrization w = 1/window means that any new observation gets the same weight as it would get in a rolling estimator with size window. Must be greater than 1.
initint, optional: Initial period. Default 10.
cutofffloat, optional: Normalized innovations exceeding this level are considered outliers. Default 2.5.

Returns:

meannp.ndarray: Exponentially weighted mean.
volnp.ndarray: Exponentially weighted volatility, i.e. a standard deviation.

cdxcore.rolling#

This Page