mfe.realized

mfe.realized

mfe.realized — Realized volatility measures for HFT data.

price_filter

price_filter(price: FloatArray, time: FloatArray, time_type: TimeType = TimeType.SECONDS, sampling_type: SamplingType = SamplingType.CALENDAR_TIME, sampling_interval: float | int | FloatArray = 300) -> tuple[FloatArray, FloatArray]

Filter raw tick prices to a regular grid.

Parameters:

Name	Type	Description	Default
price	(N,) array of log or raw prices (function is agnostic)		required
time	FloatArray		required
time_type	how timestamps are encoded		SECONDS
sampling_type	sampling scheme		CALENDAR_TIME
sampling_interval	float | int | FloatArray	CalendarTime: seconds between samples BusinessTime: number of ticks between samples CalendarUniform / BusinessUniform: number of obs in the filtered grid Fixed: (M,) array of target times	300

Returns:

Type	Description
(filtered_price, filtered_time) — both (M,) arrays

Source code in src/mfe/realized/sampling.py

def price_filter(
    price: FloatArray,
    time: FloatArray,
    time_type: TimeType = TimeType.SECONDS,
    sampling_type: SamplingType = SamplingType.CALENDAR_TIME,
    sampling_interval: float | int | FloatArray = 300,
) -> tuple[FloatArray, FloatArray]:
    """
    Filter raw tick prices to a regular grid.

    Parameters
    ----------
    price : (N,) array of log or raw prices (function is agnostic)
    time  : (N,) timestamps in units specified by time_type
    time_type : how timestamps are encoded
    sampling_type : sampling scheme
    sampling_interval :
        - CalendarTime: seconds between samples
        - BusinessTime: number of ticks between samples
        - CalendarUniform / BusinessUniform: number of obs in the filtered grid
        - Fixed: (M,) array of target times

    Returns
    -------
    (filtered_price, filtered_time) — both (M,) arrays
    """
    price = np.asarray(price, dtype=np.float64)
    time = np.asarray(time, dtype=np.float64)

    if price.shape != time.shape:
        raise ValueError(f"price and time must have the same length, got {price.shape} vs {time.shape}")

    if sampling_type == SamplingType.CALENDAR_TIME:
        return _sample_calendar_time(price, time, float(sampling_interval))
    elif sampling_type == SamplingType.BUSINESS_TIME:
        return _sample_business_time(price, time, int(sampling_interval))
    elif sampling_type == SamplingType.CALENDAR_UNIFORM:
        return _sample_calendar_uniform(price, time, int(sampling_interval))
    elif sampling_type == SamplingType.BUSINESS_UNIFORM:
        return _sample_business_uniform(price, time, int(sampling_interval))
    elif sampling_type == SamplingType.FIXED:
        target_times = np.asarray(sampling_interval, dtype=np.float64)
        return _sample_fixed(price, time, target_times)
    else:
        raise ValueError(f"Unknown sampling_type: {sampling_type}")

returns_from_prices

returns_from_prices(price: FloatArray, log: bool = True) -> FloatArray

Compute log or simple returns from a price series.

Parameters:

Name	Type	Description	Default
price	(M,) filtered price array		required
log	bool		True

Returns:

Type	Description
(M - 1,) return array

Source code in src/mfe/realized/sampling.py

def returns_from_prices(price: FloatArray, log: bool = True) -> FloatArray:
    """
    Compute log or simple returns from a price series.

    Parameters
    ----------
    price : (M,) filtered price array
    log   : if True (default), use log-price differences

    Returns
    -------
    (M - 1,) return array
    """
    price = np.asarray(price, dtype=np.float64)
    if log:
        return np.diff(np.log(price))
    else:
        return np.diff(price) / price[:-1]

refresh_time

refresh_time(prices: list[FloatArray], times: list[FloatArray]) -> tuple[list[FloatArray], FloatArray]

Synchronize K asynchronous price series via refresh-time sampling (Barndorff-Nielsen et al. 2011).

For two assets this is O(N1 + N2) and vectorized. For K > 2 this loops over assets — TODO: Cython for K > 10.

Parameters:

Name	Type	Description	Default
prices	list of K (N_k,) price arrays		required
times	list[FloatArray]		required

Returns:

Name	Type	Description
sync_prices	list of K (M,) synchronized price arrays
sync_times	(M,) refresh times

Source code in src/mfe/realized/sampling.py

def refresh_time(
    prices: list[FloatArray],
    times: list[FloatArray],
) -> tuple[list[FloatArray], FloatArray]:
    """
    Synchronize K asynchronous price series via refresh-time sampling
    (Barndorff-Nielsen et al. 2011).

    For two assets this is O(N1 + N2) and vectorized.
    For K > 2 this loops over assets — TODO: Cython for K > 10.

    Parameters
    ----------
    prices : list of K (N_k,) price arrays
    times  : list of K (N_k,) time arrays (same units)

    Returns
    -------
    sync_prices : list of K (M,) synchronized price arrays
    sync_times  : (M,) refresh times
    """
    K = len(prices)
    if K < 2:
        raise ValueError("Need at least 2 assets for refresh-time synchronization.")

    # Initial refresh time: first time all assets have a quote
    t_start = max(t[0] for t in times)

    # Build the synchronized grid iteratively
    sync_times = []
    current_idx = [np.searchsorted(times[k], t_start, side="right") - 1 for k in range(K)]
    current_idx = [max(0, i) for i in current_idx]

    while True:
        # Current refresh time = max of current "last trade" times per asset
        t_refresh = max(times[k][current_idx[k]] for k in range(K))
        sync_times.append(t_refresh)

        # Advance each asset to the first tick at or after t_refresh
        new_idx = []
        for k in range(K):
            i = np.searchsorted(times[k], t_refresh, side="left")
            new_idx.append(i)

        # Check if we've exhausted any asset
        if any(new_idx[k] >= len(times[k]) for k in range(K)):
            break

        current_idx = new_idx

    sync_times_arr = np.array(sync_times, dtype=np.float64)
    sync_prices = []
    for k in range(K):
        idx = np.searchsorted(times[k], sync_times_arr, side="right") - 1
        idx = np.clip(idx, 0, len(prices[k]) - 1)
        sync_prices.append(prices[k][idx])

    return sync_prices, sync_times_arr

realized_variance

realized_variance(returns: FloatArray, subsamples: int = 1) -> RealizedResult

Standard realized variance: sum of squared returns.

Parameters:

Name	Type	Description	Default
returns	FloatArray		required
subsamples	number of sub-grids for sub-sampling bias correction		1

Returns:

Type	Description
RealizedResult with .value = RV and .subsampled_value = sub-sampled RV

Source code in src/mfe/realized/variance.py

def realized_variance(
    returns: FloatArray,
    subsamples: int = 1,
) -> RealizedResult:
    """
    Standard realized variance: sum of squared returns.

    Parameters
    ----------
    returns    : (M,) log-return array
    subsamples : number of sub-grids for sub-sampling bias correction

    Returns
    -------
    RealizedResult with .value = RV and .subsampled_value = sub-sampled RV
    """
    r = np.asarray(returns, dtype=np.float64)
    rv = float(np.sum(r ** 2))

    rv_ss = None
    if subsamples > 1:
        ss_rvs = []
        for s in range(subsamples):
            r_sub = r[s::subsamples]
            ss_rvs.append(float(np.sum(r_sub ** 2)))
        rv_ss = float(np.mean(ss_rvs)) * subsamples  # scale back to full-sample

    return RealizedResult(
        value=rv,
        subsampled_value=rv_ss,
        n_returns=len(r),
    )

realized_bipower_variation

realized_bipower_variation(returns: FloatArray, skip: int = 0, subsamples: int = 1) -> RealizedResult

Realized bipower variation (BPV) with optional skip-k extension.

BPV = mu_1^{-2} * sum_{t=skip+2}^{T} |r_t| * |r_{t-skip-1}|

Parameters:

Name	Type	Description	Default
returns	FloatArray		required
skip	int		0
subsamples	sub-sampling replications for bias correction		1

Returns:

Type	Description
RealizedResult	.value = BPV .debiased_value = BPV * m/(m - skip - 1) where m = number of returns used

Source code in src/mfe/realized/variance.py

def realized_bipower_variation(
    returns: FloatArray,
    skip: int = 0,
    subsamples: int = 1,
) -> RealizedResult:
    """
    Realized bipower variation (BPV) with optional skip-k extension.

    BPV = mu_1^{-2} * sum_{t=skip+2}^{T} |r_t| * |r_{t-skip-1}|

    Parameters
    ----------
    returns    : (M,) log-return array
    skip       : number of returns to skip between the two absolute returns
    subsamples : sub-sampling replications for bias correction

    Returns
    -------
    RealizedResult
        .value           = BPV
        .debiased_value  = BPV * m/(m - skip - 1) where m = number of returns used
    """
    r = np.asarray(returns, dtype=np.float64)
    bpv, m = _bpv_core(r, skip)

    debiased = bpv * m / (m - skip - 1) if m > skip + 1 else np.nan

    bpv_ss = None
    if subsamples > 1:
        ss_vals = []
        for s in range(subsamples):
            r_sub = r[s::subsamples]
            val, m_sub = _bpv_core(r_sub, skip)
            ss_vals.append(val * subsamples)
        bpv_ss = float(np.mean(ss_vals))

    return RealizedResult(
        value=bpv,
        subsampled_value=bpv_ss,
        debiased_value=float(debiased),
        n_returns=len(r),
    )

realized_med_variance

realized_med_variance(returns: FloatArray) -> RealizedResult

Median realized variance: robust to jumps.

MedRV = (pi / (6 - 4*sqrt(3) + pi)) * (M/(M-2)) * sum_{t=2}^{T-1} median(|r_{t-1}|, |r_t|, |r_{t+1}|)^2

Vectorized: uses np.partition (O(N), not O(N log N)) to find the median of each triplet without sorting. ~3x faster than the column_stack approach.

Source code in src/mfe/realized/variance.py

def realized_med_variance(returns: FloatArray) -> RealizedResult:
    """
    Median realized variance: robust to jumps.

    MedRV = (pi / (6 - 4*sqrt(3) + pi)) * (M/(M-2)) *
            sum_{t=2}^{T-1} median(|r_{t-1}|, |r_t|, |r_{t+1}|)^2

    Vectorized: uses np.partition (O(N), not O(N log N)) to find the median
    of each triplet without sorting. ~3x faster than the column_stack approach.
    """
    r = np.asarray(returns, dtype=np.float64)
    M = len(r)
    absr = np.abs(r)

    a0 = absr[:-2]
    a1 = absr[1:-1]
    a2 = absr[2:]

    if _HAS_CYTHON:
        raw_sum = float(_medvar_cy(np.ascontiguousarray(absr, dtype=np.float64)))
    else:
        triplets = np.stack([a0, a1, a2], axis=1)
        partitioned = np.partition(triplets, kth=1, axis=1)
        raw_sum = float(np.sum(partitioned[:, 1] ** 2))

    pi = np.pi
    scale = (pi / (6 - 4 * np.sqrt(3) + pi)) * (M / (M - 2))
    med_rv = float(scale * raw_sum)

    return RealizedResult(value=med_rv, n_returns=M)

realized_min_variance

realized_min_variance(returns: FloatArray) -> RealizedResult

Min realized variance: minimum of adjacent pairs of squared returns.

MinRV = (pi / (pi - 2)) * (M/(M-1)) * sum_{t=1}^{T-1} min(|r_t|, |r_{t+1}|)^2

Source code in src/mfe/realized/variance.py

def realized_min_variance(returns: FloatArray) -> RealizedResult:
    """
    Min realized variance: minimum of adjacent pairs of squared returns.

    MinRV = (pi / (pi - 2)) * (M/(M-1)) * sum_{t=1}^{T-1} min(|r_t|, |r_{t+1}|)^2
    """
    r = np.asarray(returns, dtype=np.float64)
    M = len(r)
    absr = np.abs(r)

    pairs = np.column_stack([absr[:-1], absr[1:]])
    min_sq = np.min(pairs, axis=1) ** 2
    pi = np.pi
    scale = (pi / (pi - 2)) * (M / (M - 1))
    min_rv = float(scale * np.sum(min_sq))

    return RealizedResult(value=min_rv, n_returns=M)

realized_preaveraged_variance

realized_preaveraged_variance(returns: FloatArray, theta: float = 0.8) -> RealizedResult

Pre-averaged realized variance (Jacod et al. 2009).

Uses a linear pre-averaging kernel g(x) = min(x, 1-x) with block size k_n = floor(theta * sqrt(n)).

This estimator is consistent even under microstructure noise.

Parameters:

Name	Type	Description	Default
returns	(M,) log-return array		required
theta	float		0.8

Source code in src/mfe/realized/variance.py

def realized_preaveraged_variance(
    returns: FloatArray,
    theta: float = 0.8,
) -> RealizedResult:
    """
    Pre-averaged realized variance (Jacod et al. 2009).

    Uses a linear pre-averaging kernel g(x) = min(x, 1-x) with block size
    k_n = floor(theta * sqrt(n)).

    This estimator is consistent even under microstructure noise.

    Parameters
    ----------
    returns : (M,) log-return array
    theta   : tuning parameter controlling block size (default 0.8)
    """
    r = np.asarray(returns, dtype=np.float64)
    n = len(r)
    kn = max(2, int(np.floor(theta * np.sqrt(n))))

    # g(x) = min(x, 1-x) evaluated at x = i/kn for i=1..kn-1
    i_vals = np.arange(1, kn)
    g = np.minimum(i_vals / kn, 1 - i_vals / kn)
    g_sq_sum = float(np.sum(g ** 2))
    psi2 = g_sq_sum / kn  # psi_2 normalization constant

    # Pre-average
    pre_avg = np.zeros(n - kn + 1, dtype=np.float64)
    for j in range(kn - 1):
        pre_avg[: n - kn + 1] += g[j] * r[j : n - kn + 1 + j]

    pv = float(np.sum(pre_avg ** 2)) / (kn * psi2)

    # Noise bias correction (uses realized variance at fine scale)
    rv_fine = float(np.sum(r ** 2))
    bias = (kn / 2) * psi2 * rv_fine
    pv_corrected = pv - bias / kn

    return RealizedResult(
        value=pv_corrected,
        n_returns=n,
        diagnostics={"kn": kn, "theta": theta, "psi2": psi2},
    )

realized_semivariance

realized_semivariance(returns: FloatArray) -> tuple[RealizedResult, RealizedResult]

Decompose RV into positive and negative semivariance.

RS+ = sum_{r > 0} r^2, RS- = sum_{r < 0} r^2

Returns (rs_pos, rs_neg).

Source code in src/mfe/realized/variance.py

def realized_semivariance(
    returns: FloatArray,
) -> tuple[RealizedResult, RealizedResult]:
    """
    Decompose RV into positive and negative semivariance.

    RS+ = sum_{r > 0} r^2,  RS- = sum_{r < 0} r^2

    Returns (rs_pos, rs_neg).
    """
    r = np.asarray(returns, dtype=np.float64)
    rs_pos = float(np.sum(r[r > 0] ** 2))
    rs_neg = float(np.sum(r[r < 0] ** 2))
    return (
        RealizedResult(value=rs_pos, n_returns=int(np.sum(r > 0))),
        RealizedResult(value=rs_neg, n_returns=int(np.sum(r < 0))),
    )

realized_kernel

realized_kernel(returns: FloatArray, kernel_type: KernelType = KernelType.PARZEN, bandwidth: int | None = None, jitter: bool = True) -> RealizedKernelResult

Realized kernel estimator for quadratic variation.

Parameters:

Name	Type	Description	Default
returns	FloatArray		required
kernel_type	KernelType		PARZEN
bandwidth	int | None		None
jitter	bool	as in the BNHLS paper; adds a small fraction of the noise variance to handle the boundary bias	True

Returns:

Type	Description
RealizedKernelResult

Source code in src/mfe/realized/kernel.py

def realized_kernel(
    returns: FloatArray,
    kernel_type: KernelType = KernelType.PARZEN,
    bandwidth: int | None = None,
    jitter: bool = True,
) -> RealizedKernelResult:
    """
    Realized kernel estimator for quadratic variation.

    Parameters
    ----------
    returns      : (M,) log-return array (already filtered/sampled)
    kernel_type  : which kernel weight function to use
    bandwidth    : H; if None, uses automatic selector
    jitter       : if True, apply end-point jittering (noise correction)
                   as in the BNHLS paper; adds a small fraction of the
                   noise variance to handle the boundary bias

    Returns
    -------
    RealizedKernelResult
    """
    r = np.asarray(returns, dtype=np.float64)
    n = len(r)

    # Noise variance (needed for bandwidth selection and jitter)
    noise_var = estimate_noise_variance(r)

    if bandwidth is None:
        H = select_bandwidth(r, kernel_type=kernel_type, noise_variance=noise_var)
    else:
        H = int(bandwidth)

    if H >= n:
        warnings.warn(
            f"Bandwidth H={H} >= n={n}; clipping to n//2.",
            RuntimeWarning,
            stacklevel=2,
        )
        H = n // 2

    # Compute autocovariances
    if _HAS_CYTHON:
        gamma = _acov_fast(r, H)
    else:
        gamma = _autocovariance_numpy(r, H)

    # Kernel weights
    w = _kernel_weights(kernel_type, H)

    # RK = gamma_0 + 2 * sum_{h=1}^{H} k(h/(H+1)) * gamma_h
    rk = gamma[0] + 2.0 * float(w[1:] @ gamma[1:])

    # End-point (jitter) correction: adjusts for noise at boundaries
    # See BNHLS eq. (29); adds 2 * noise_var
    rk_adjusted = rk - 2.0 * noise_var if jitter else rk

    rk_adjusted = max(rk_adjusted, 0.0)  # enforce positivity

    return RealizedKernelResult(
        rk=rk,
        rk_adjusted=rk_adjusted,
        bandwidth=H,
        noise_variance=noise_var,
        iq_lower_bound=0.0,  # populated by caller if needed
        kernel_type=kernel_type,
        n_returns=n,
    )

select_bandwidth

select_bandwidth(returns: FloatArray, kernel_type: KernelType = KernelType.PARZEN, noise_variance: float | None = None, iq_lower_bound: float | None = None) -> int

Optimal bandwidth H for the realized kernel.

H* = c_star * xi^{4/5} * n^{3/5}

where xi = noise_variance / sqrt(IQ), and c_star depends on the kernel.

Source code in src/mfe/realized/kernel.py

def select_bandwidth(
    returns: FloatArray,
    kernel_type: KernelType = KernelType.PARZEN,
    noise_variance: float | None = None,
    iq_lower_bound: float | None = None,
) -> int:
    """
    Optimal bandwidth H for the realized kernel.

    H* = c_star * xi^{4/5} * n^{3/5}

    where xi = noise_variance / sqrt(IQ), and c_star depends on the kernel.
    """
    r = np.asarray(returns, dtype=np.float64)
    n = len(r)

    if noise_variance is None:
        noise_variance = estimate_noise_variance(r)

    # Lower bound for IQ: use tripower variation as proxy
    if iq_lower_bound is None:
        from mfe.realized.quarticity import realized_quarticity
        iq_lower_bound = realized_quarticity(r).value

    # c_star depends on kernel (from Table 1 of BNHLS 2009)
    c_star_map = {
        KernelType.PARZEN: 3.51,
        KernelType.BARTLETT: 2.16,
        KernelType.TUKEY_HANNING: 3.68,
        KernelType.CUBIC: 3.71,
        KernelType.EPANECHNIKOV: 3.28,
        KernelType.FLAT_TOP: 2.78,
    }
    c_star = c_star_map.get(kernel_type, 3.51)

    xi = noise_variance / max(iq_lower_bound ** 0.5, 1e-30)
    H = max(1, int(np.round(c_star * (xi ** 0.4) * (n ** 0.6))))

    return H

realized_quarticity

realized_quarticity(returns: FloatArray) -> RealizedResult

Realized quarticity: (n/3) * sum r_t^4

Consistent estimator of integrated quarticity IQ = int_0^1 sigma_t^4 dt.

Source code in src/mfe/realized/quarticity.py

def realized_quarticity(returns: FloatArray) -> RealizedResult:
    """
    Realized quarticity: (n/3) * sum r_t^4

    Consistent estimator of integrated quarticity IQ = int_0^1 sigma_t^4 dt.
    """
    r = np.asarray(returns, dtype=np.float64)
    n = len(r)
    rq = float((n / 3) * np.sum(r ** 4))
    return RealizedResult(value=rq, n_returns=n)

realized_tripower_quarticity

realized_tripower_quarticity(returns: FloatArray) -> RealizedResult

Tripower quarticity — robust to occasional jumps.

TPQ = n * mu_{4/3}^{-3} * mean(|r_{t-2}|^{4/3} |r_{t-1}|^{4/3} |r_t|^{4/3})

Source code in src/mfe/realized/quarticity.py

def realized_tripower_quarticity(returns: FloatArray) -> RealizedResult:
    """
    Tripower quarticity — robust to occasional jumps.

    TPQ = n * mu_{4/3}^{-3} * mean(|r_{t-2}|^{4/3} |r_{t-1}|^{4/3} |r_t|^{4/3})
    """
    r = np.asarray(returns, dtype=np.float64)
    n = len(r)
    absr = np.abs(r)
    tpq = float(
        n * (_MU43 ** -3) * np.mean(absr[:-2] ** (4 / 3) * absr[1:-1] ** (4 / 3) * absr[2:] ** (4 / 3))
    )
    return RealizedResult(value=tpq, n_returns=n)

bns_jump_test

bns_jump_test(returns: FloatArray, alpha: float = 0.05) -> JumpTestResult

Barndorff-Nielsen & Shephard (2006) jump test based on the ratio RV/BPV.

Z = sqrt(n) * (RV/BPV - 1) / sqrt(omega_hat)

Under the null of no jumps, Z -> N(0, 1).

Parameters:

Name	Type	Description	Default
returns	(M,) log-return array		required
alpha	float		0.05

Returns:

Type	Description
JumpTestResult

Source code in src/mfe/realized/jumps.py

def bns_jump_test(
    returns: FloatArray,
    alpha: float = 0.05,
) -> JumpTestResult:
    """
    Barndorff-Nielsen & Shephard (2006) jump test based on the ratio RV/BPV.

    Z = sqrt(n) * (RV/BPV - 1) / sqrt(omega_hat)

    Under the null of no jumps, Z -> N(0, 1).

    Parameters
    ----------
    returns : (M,) log-return array
    alpha   : significance level

    Returns
    -------
    JumpTestResult
    """
    r = np.asarray(returns, dtype=np.float64)
    n = len(r)

    rv = realized_variance(r).value
    bpv = realized_bipower_variation(r).value
    tpq = realized_tripower_quarticity(r).value

    # Consistent estimate of asymptotic variance (BN-S 2006, Theorem 1)
    # omega = (pi^2/4 + pi - 5) * max(1, TQ/BPV^2)
    pi = np.pi
    omega_hat = (pi ** 2 / 4 + pi - 5) * max(1.0, tpq / max(bpv ** 2, 1e-300))

    z_stat = float(np.sqrt(n) * (rv / max(bpv, 1e-300) - 1) / np.sqrt(max(omega_hat, 1e-300)))
    p_val = float(2 * (1 - stats.norm.cdf(abs(z_stat))))

    jump_var = max(rv - bpv, 0.0)
    return JumpTestResult(
        statistic=z_stat,
        p_value=p_val,
        jump_variation=jump_var,
        continuous_variation=bpv,
        total_variation=rv,
        significant=(p_val < alpha),
    )

estimate_noise_variance

estimate_noise_variance(returns: FloatArray, method: str = 'bandi-russell') -> float

Estimate the microstructure noise variance omega^2.

Parameters:

Name	Type	Description	Default
returns	(M,) log-return array at the finest available frequency		required
method	str		'bandi-russell'

Returns:

Type	Description
float — noise variance estimate (>= 0)

Source code in src/mfe/realized/noise.py

def estimate_noise_variance(
    returns: FloatArray,
    method: str = "bandi-russell",
) -> float:
    """
    Estimate the microstructure noise variance omega^2.

    Parameters
    ----------
    returns : (M,) log-return array at the finest available frequency
    method  : "bandi-russell" (default) or "zma"

    Returns
    -------
    float — noise variance estimate (>= 0)
    """
    r = np.asarray(returns, dtype=np.float64)

    if method == "bandi-russell":
        return _noise_bandi_russell(r)
    elif method == "zma":
        return _noise_zma(r)
    else:
        raise ValueError(f"Unknown noise estimation method: {method}")

realized_covariance

realized_covariance(returns: FloatArray) -> RealizedCovarianceResult

Standard realized covariance matrix from synchronous returns.

Parameters:

Name	Type	Description	Default
returns	(T, K) matrix of synchronous log-returns		required

Returns:

Type	Description
RealizedCovarianceResult with .cov = (K, K) realized covariance matrix

Source code in src/mfe/realized/covariance.py

def realized_covariance(
    returns: FloatArray,
) -> RealizedCovarianceResult:
    """
    Standard realized covariance matrix from synchronous returns.

    Parameters
    ----------
    returns : (T, K) matrix of synchronous log-returns

    Returns
    -------
    RealizedCovarianceResult with .cov = (K, K) realized covariance matrix
    """
    r = np.asarray(returns, dtype=np.float64)
    if r.ndim == 1:
        r = r[:, None]
    T, K = r.shape

    cov = r.T @ r  # (K, K) — NOT divided by T, this is the quadratic variation

    return RealizedCovarianceResult(
        cov=cov,
        method="synchronous",
        n_assets=K,
        n_returns=T,
    )

realized_correlation

realized_correlation(returns: FloatArray) -> FloatArray

Realized correlation matrix from synchronous returns.

Returns (K, K) correlation matrix.

Source code in src/mfe/realized/covariance.py

def realized_correlation(returns: FloatArray) -> FloatArray:
    """
    Realized correlation matrix from synchronous returns.

    Returns (K, K) correlation matrix.
    """
    res = realized_covariance(returns)
    cov = res.cov
    d = np.sqrt(np.diag(cov))
    d_inv = np.where(d > 0, 1.0 / d, 0.0)
    return d_inv[:, None] * cov * d_inv[None, :]

realized_hayashi_yoshida

realized_hayashi_yoshida(prices: list[FloatArray], times: list[FloatArray]) -> RealizedCovarianceResult

Hayashi-Yoshida realized covariance for K non-synchronously observed assets.

Hayashi, T. & Yoshida, N. (2005): "On Covariance Estimation of Non-Synchronously Observed Diffusion Processes", Bernoulli.

Parameters:

Name	Type	Description	Default
prices	list of K price arrays (lengths can differ)		required
times	list[FloatArray]		required

Returns:

Type	Description
RealizedCovarianceResult with (K, K) covariance matrix.	.method = "hayashi-yoshida"

Notes

Diagonal elements are the standard realized variance of each asset (computed from their own tick data, so no synchronization needed).

K > 2 assets: implemented as O(K^2) bivariate calls. The MATLAB mfe-toolbox has a TODO here for the general case — we implement it.

Source code in src/mfe/realized/covariance.py

def realized_hayashi_yoshida(
    prices: list[FloatArray],
    times: list[FloatArray],
) -> RealizedCovarianceResult:
    """
    Hayashi-Yoshida realized covariance for K non-synchronously observed assets.

    Hayashi, T. & Yoshida, N. (2005): "On Covariance Estimation of
    Non-Synchronously Observed Diffusion Processes", Bernoulli.

    Parameters
    ----------
    prices : list of K price arrays (lengths can differ)
    times  : list of K timestamp arrays

    Returns
    -------
    RealizedCovarianceResult with (K, K) covariance matrix.
        .method = "hayashi-yoshida"

    Notes
    -----
    Diagonal elements are the standard realized variance of each asset
    (computed from their own tick data, so no synchronization needed).

    K > 2 assets: implemented as O(K^2) bivariate calls. The MATLAB
    mfe-toolbox has a TODO here for the general case — we implement it.
    """
    K = len(prices)
    if K < 2:
        raise ValueError("Need at least 2 assets.")
    if len(times) != K:
        raise ValueError(f"len(times)={len(times)} != len(prices)={K}")

    cov = np.zeros((K, K), dtype=np.float64)

    # Diagonal: realized variance from each asset's own tick data
    for k in range(K):
        r_k = np.diff(np.log(np.asarray(prices[k], dtype=np.float64)))
        cov[k, k] = float(np.sum(r_k ** 2))

    # Off-diagonal: HY estimator for each pair
    for i in range(K):
        for j in range(i + 1, K):
            hy = _hy_bivariate(
                np.asarray(prices[i], dtype=np.float64),
                np.asarray(times[i], dtype=np.float64),
                np.asarray(prices[j], dtype=np.float64),
                np.asarray(times[j], dtype=np.float64),
            )
            cov[i, j] = hy
            cov[j, i] = hy

    n_returns = min(len(p) - 1 for p in prices)

    return RealizedCovarianceResult(
        cov=cov,
        method="hayashi-yoshida",
        n_assets=K,
        n_returns=n_returns,
    )

realized_covariance_refresh_time

realized_covariance_refresh_time(prices: list[FloatArray], times: list[FloatArray]) -> RealizedCovarianceResult

Realized covariance using refresh-time synchronization.

Barndorff-Nielsen, Hansen, Lunde & Shephard (2011).

Synchronizes K asynchronous tick streams via refresh time, then computes the standard realized covariance on the synchronized returns.

Less efficient than HY (more data loss from synchronization) but gives a positive semi-definite matrix by construction.

Source code in src/mfe/realized/covariance.py

def realized_covariance_refresh_time(
    prices: list[FloatArray],
    times: list[FloatArray],
) -> RealizedCovarianceResult:
    """
    Realized covariance using refresh-time synchronization.

    Barndorff-Nielsen, Hansen, Lunde & Shephard (2011).

    Synchronizes K asynchronous tick streams via refresh time, then computes
    the standard realized covariance on the synchronized returns.

    Less efficient than HY (more data loss from synchronization) but gives a
    positive semi-definite matrix by construction.
    """
    sync_prices, sync_times = refresh_time(prices, times)
    sync_returns = np.column_stack([
        np.diff(np.log(p)) for p in sync_prices
    ])
    return realized_covariance(sync_returns)

realized_range

realized_range(high: FloatArray, low: FloatArray, open_: FloatArray | None = None, close: FloatArray | None = None) -> RealizedRangeResult

Realized range estimator from sub-interval high/low prices.

Parameters:

Name	Type	Description	Default
high	FloatArray		required
low	FloatArray		required
open_	(M,) optional — open log-price (used for Yang-Zhang correction)		None
close	(M,) optional — close log-price (used for Yang-Zhang correction)		None

Returns:

Type	Description
RealizedRangeResult	.value — realized range estimate of IV .n_intervals — M .efficiency — ~1.67 vs RV under Brownian motion

Notes

Input prices can be raw or log — the function uses log(high) - log(low) which equals log(high/low) regardless of whether inputs are already logs. If inputs are already log-prices, pass them directly; the formula is the same either way (log(e^h) - log(e^l) = h - l).

Source code in src/mfe/realized/range_.py

def realized_range(
    high: FloatArray,
    low: FloatArray,
    open_: FloatArray | None = None,
    close: FloatArray | None = None,
) -> RealizedRangeResult:
    """
    Realized range estimator from sub-interval high/low prices.

    Parameters
    ----------
    high  : (M,) highest log-price in each sub-interval (or raw price)
    low   : (M,) lowest log-price in each sub-interval
    open_ : (M,) optional — open log-price (used for Yang-Zhang correction)
    close : (M,) optional — close log-price (used for Yang-Zhang correction)

    Returns
    -------
    RealizedRangeResult
        .value — realized range estimate of IV
        .n_intervals — M
        .efficiency — ~1.67 vs RV under Brownian motion

    Notes
    -----
    Input prices can be raw or log — the function uses log(high) - log(low)
    which equals log(high/low) regardless of whether inputs are already logs.
    If inputs are already log-prices, pass them directly; the formula is
    the same either way (log(e^h) - log(e^l) = h - l).
    """
    H = np.asarray(high, dtype=np.float64)
    L = np.asarray(low, dtype=np.float64)

    if H.shape != L.shape:
        raise ValueError(f"high and low must have the same shape, got {H.shape} vs {L.shape}")

    if np.any(H < L):
        raise ValueError("high must be >= low in every interval")

    M = len(H)
    log_range_sq = (np.log(H) - np.log(L)) ** 2

    rr = float(np.sum(log_range_sq)) / _4LOG2

    # Theoretical efficiency of range vs squared-return: ~1.67 under BM
    efficiency = 1.6704  # exact: 1 / (4 log 2) * pi^2/2 * (1 - 2/pi)^{-1}; approx

    return RealizedRangeResult(value=rr, n_intervals=M, efficiency=efficiency)

realized_range_from_ticks

realized_range_from_ticks(price: FloatArray, time: FloatArray, interval_seconds: float = 300.0) -> RealizedRangeResult

Compute realized range from raw tick data by aggregating into OHLC bars.

Parameters:

Name	Type	Description	Default
price	FloatArray		required
time	FloatArray		required
interval_seconds	bar width in seconds (default 5 minutes)		300.0

Source code in src/mfe/realized/range_.py

def realized_range_from_ticks(
    price: FloatArray,
    time: FloatArray,
    interval_seconds: float = 300.0,
) -> RealizedRangeResult:
    """
    Compute realized range from raw tick data by aggregating into OHLC bars.

    Parameters
    ----------
    price            : (N,) raw tick prices
    time             : (N,) timestamps in seconds
    interval_seconds : bar width in seconds (default 5 minutes)
    """
    price = np.asarray(price, dtype=np.float64)
    time  = np.asarray(time, dtype=np.float64)

    t_start = time[0]
    t_end   = time[-1]

    # Build bar boundaries
    edges = np.arange(t_start, t_end + interval_seconds, interval_seconds)
    M = len(edges) - 1

    highs  = np.empty(M, dtype=np.float64)
    lows   = np.empty(M, dtype=np.float64)
    valid  = np.zeros(M, dtype=bool)

    for j in range(M):
        mask = (time >= edges[j]) & (time < edges[j + 1])
        if np.any(mask):
            bar_prices = price[mask]
            highs[j] = np.max(bar_prices)
            lows[j]  = np.min(bar_prices)
            valid[j] = True

    return realized_range(highs[valid], lows[valid])

msrv

msrv(returns: FloatArray, n_scales: int | None = None) -> MSRVResult

Multi-Scale Realized Variance (MSRV) — Zhang (2006).

Combines J sub-sampled RVs with optimally chosen weights to achieve the best rate of convergence under i.i.d. microstructure noise.

Parameters:

Name	Type	Description	Default
returns	FloatArray		required
n_scales	number of scales J; if None uses min(N^{1/2}, 30)		None

Source code in src/mfe/realized/tsrv.py

def msrv(
    returns: FloatArray,
    n_scales: int | None = None,
) -> MSRVResult:
    """
    Multi-Scale Realized Variance (MSRV) — Zhang (2006).

    Combines J sub-sampled RVs with optimally chosen weights to achieve
    the best rate of convergence under i.i.d. microstructure noise.

    Parameters
    ----------
    returns  : (N,) all-tick log-returns
    n_scales : number of scales J; if None uses min(N^{1/2}, 30)
    """
    r = np.asarray(returns, dtype=np.float64)
    N = len(r)

    J = n_scales if n_scales is not None else min(30, max(2, int(np.sqrt(N))))

    # Sub-sampled RVs at scales K = 1, 2, ..., J
    rv_scales = np.empty(J, dtype=np.float64)
    n_returns = np.empty(J, dtype=np.float64)

    for j, K in enumerate(range(1, J + 1)):
        rv_grids = [float(np.sum(r[s::K] ** 2)) for s in range(K)]
        rv_scales[j] = float(np.mean(rv_grids))
        n_returns[j] = N / K

    # Optimal weights from Zhang (2006) eq. (3.13)
    # w_K = 12 * K * (K - 1/2) * (1/n_K) / (J(J+1)(2J+1))
    K_arr = np.arange(1, J + 1, dtype=float)
    denom = J * (J + 1) * (2 * J + 1)
    weights = 12 * K_arr * (K_arr - 0.5) / (n_returns * denom)

    msrv_val = float(weights @ rv_scales)
    rv_fast = float(np.sum(r ** 2))

    return MSRVResult(
        msrv=msrv_val,
        rv_fast=rv_fast,
        n_scales=J,
        weights=weights,
        rv_per_scale=rv_scales,
    )

realized_quantile_variance

realized_quantile_variance(returns: FloatArray, tau: float = 0.5) -> RealizedQuantileVarResult

Realized quantile variance — jump-robust quadratic variation estimator.

Parameters:

Name	Type	Description	Default
returns	(T,) log-return array		required
tau	float	Lower tau → more jump-robust, less efficient. tau = 1 → equivalent to realized variance (no truncation).	0.5

Returns:

Type	Description
RealizedQuantileVarResult	.value — RQV estimate of integrated variance

Source code in src/mfe/realized/quantile_var.py

def realized_quantile_variance(
    returns: FloatArray,
    tau: float = 0.50,
) -> RealizedQuantileVarResult:
    """
    Realized quantile variance — jump-robust quadratic variation estimator.

    Parameters
    ----------
    returns : (T,) log-return array
    tau     : quantile probability in (0, 1); default 0.50 (median)
              Lower tau → more jump-robust, less efficient.
              tau = 1 → equivalent to realized variance (no truncation).

    Returns
    -------
    RealizedQuantileVarResult
        .value — RQV estimate of integrated variance
    """
    if not 0 < tau < 1:
        raise ValueError(f"tau must be in (0, 1), got {tau}")

    r = np.asarray(returns, dtype=np.float64)
    T = len(r)
    r2 = r ** 2

    # Quantile cutoff on squared returns
    q_tau = float(np.quantile(r2, tau))

    # Truncated sum: only squared returns below q_tau
    mask = r2 <= q_tau
    n_used = int(np.sum(mask))
    truncated_sum = float(np.sum(r2[mask]))

    if n_used == 0:
        return RealizedQuantileVarResult(
            value=np.nan, tau=tau, n_returns=T, n_truncated=0, quantile_cutoff=q_tau
        )

    # Raw estimate: scale so that the mean over ALL T observations is calibrated
    raw = truncated_sum / T

    # Calibration constant
    c = _calibration_constant(tau)

    rqv = raw / c if c > 0 else raw

    return RealizedQuantileVarResult(
        value=rqv,
        tau=tau,
        n_returns=T,
        n_truncated=n_used,
        quantile_cutoff=q_tau,
    )

realized_multivariate_kernel

realized_multivariate_kernel(returns: FloatArray, kernel_type: KernelType = KernelType.PARZEN, bandwidth: int | None = None, jitter: bool = True) -> MultivariateKernelResult

Multivariate realized kernel for the full (K, K) covariance matrix.

Parameters:

Name	Type	Description	Default
returns	FloatArray		required
kernel_type	kernel weight function		PARZEN
bandwidth	int | None		None
jitter	bool	2 * diag(noise_var) from the diagonal)	True

Returns:

Type	Description
MultivariateKernelResult	.rk — raw (K, K) realized kernel matrix .rk_adjusted — noise-corrected version (PSD enforced)

Source code in src/mfe/realized/multivariate_kernel.py

def realized_multivariate_kernel(
    returns: FloatArray,
    kernel_type: KernelType = KernelType.PARZEN,
    bandwidth: int | None = None,
    jitter: bool = True,
) -> MultivariateKernelResult:
    """
    Multivariate realized kernel for the full (K, K) covariance matrix.

    Parameters
    ----------
    returns     : (T, K) synchronized log-return matrix
    kernel_type : kernel weight function
    bandwidth   : H; if None, uses average of per-asset optimal bandwidths
    jitter      : if True, apply end-point noise correction (subtract
                  2 * diag(noise_var) from the diagonal)

    Returns
    -------
    MultivariateKernelResult
        .rk          — raw (K, K) realized kernel matrix
        .rk_adjusted — noise-corrected version (PSD enforced)
    """
    R = np.asarray(returns, dtype=np.float64)
    if R.ndim == 1:
        raise ValueError("returns must be (T, K) with K >= 2 for multivariate kernel")
    T, K = R.shape

    # Per-asset noise variances and optimal bandwidths
    noise_vars = np.array([estimate_noise_variance(R[:, k]) for k in range(K)])

    if bandwidth is None:
        H_per_asset = [select_bandwidth(R[:, k], kernel_type=kernel_type,
                                        noise_variance=noise_vars[k]) for k in range(K)]
        H = max(1, int(np.round(np.mean(H_per_asset))))
    else:
        H = int(bandwidth)

    if H >= T:
        warnings.warn(f"Bandwidth H={H} >= T={T}; clipping to T//2.", RuntimeWarning, stacklevel=2)
        H = T // 2

    # Kernel weights
    w = _kernel_weights(kernel_type, H)   # (H+1,)

    # Cross-autocovariance matrices Gamma_h = R[h:].T @ R[:T-h]  for h=0..H
    # Stack: (H+1, K, K) then weight-sum
    RK = np.zeros((K, K), dtype=np.float64)

    # h = 0: R.T @ R
    RK += w[0] * (R.T @ R)

    # h = 1..H: w[h] * (Gamma_h + Gamma_h.T) — symmetric kernel
    for h in range(1, H + 1):
        Gamma_h = R[h:].T @ R[:T - h]   # (K, K)
        RK += w[h] * (Gamma_h + Gamma_h.T)

    # Jitter correction: subtract 2 * diag(noise_var) per asset
    if jitter:
        noise_correction = np.diag(2.0 * noise_vars)
        rk_adj = RK - noise_correction
    else:
        rk_adj = RK.copy()

    # Enforce PSD: project onto cone of PSD matrices
    eigvals, eigvecs = np.linalg.eigh(rk_adj)
    if np.any(eigvals < 0):
        eigvals_clipped = np.maximum(eigvals, 0.0)
        rk_adj = eigvecs @ np.diag(eigvals_clipped) @ eigvecs.T

    return MultivariateKernelResult(
        rk=RK,
        rk_adjusted=rk_adj,
        bandwidth=H,
        noise_variances=noise_vars,
        kernel_type=kernel_type,
        n_returns=T,
        n_vars=K,
    )

covariance

Realized covariance estimators for multivariate HFT data.

Andersen, Bollerslev, Diebold & Labys (2003): synchronous realized covariance. Hayashi & Yoshida (2005): non-synchronous covariance estimator. Barndorff-Nielsen et al. (2011): multivariate realized kernel.

Note on Hayashi-Yoshida for K > 2: The MATLAB realized_hayashi_yoshida.m has a TODO comment for K > 2 assets. We implement the general K-asset case by applying the bivariate HY estimator to each (i, j) pair and assembling the full matrix. This is O(K^2 * max(N_i, N_j)) and correct, but not the most efficient possible implementation for large K.

realized_covariance

realized_covariance(returns: FloatArray) -> RealizedCovarianceResult

Standard realized covariance matrix from synchronous returns.

Parameters:

Name	Type	Description	Default
returns	(T, K) matrix of synchronous log-returns		required

Returns:

Type	Description
RealizedCovarianceResult with .cov = (K, K) realized covariance matrix

Source code in src/mfe/realized/covariance.py

def realized_covariance(
    returns: FloatArray,
) -> RealizedCovarianceResult:
    """
    Standard realized covariance matrix from synchronous returns.

    Parameters
    ----------
    returns : (T, K) matrix of synchronous log-returns

    Returns
    -------
    RealizedCovarianceResult with .cov = (K, K) realized covariance matrix
    """
    r = np.asarray(returns, dtype=np.float64)
    if r.ndim == 1:
        r = r[:, None]
    T, K = r.shape

    cov = r.T @ r  # (K, K) — NOT divided by T, this is the quadratic variation

    return RealizedCovarianceResult(
        cov=cov,
        method="synchronous",
        n_assets=K,
        n_returns=T,
    )

realized_correlation

realized_correlation(returns: FloatArray) -> FloatArray

Realized correlation matrix from synchronous returns.

Returns (K, K) correlation matrix.

Source code in src/mfe/realized/covariance.py

def realized_correlation(returns: FloatArray) -> FloatArray:
    """
    Realized correlation matrix from synchronous returns.

    Returns (K, K) correlation matrix.
    """
    res = realized_covariance(returns)
    cov = res.cov
    d = np.sqrt(np.diag(cov))
    d_inv = np.where(d > 0, 1.0 / d, 0.0)
    return d_inv[:, None] * cov * d_inv[None, :]

realized_hayashi_yoshida

realized_hayashi_yoshida(prices: list[FloatArray], times: list[FloatArray]) -> RealizedCovarianceResult

Hayashi-Yoshida realized covariance for K non-synchronously observed assets.

Hayashi, T. & Yoshida, N. (2005): "On Covariance Estimation of Non-Synchronously Observed Diffusion Processes", Bernoulli.

Parameters:

Name	Type	Description	Default
prices	list of K price arrays (lengths can differ)		required
times	list[FloatArray]		required

Returns:

Type	Description
RealizedCovarianceResult with (K, K) covariance matrix.	.method = "hayashi-yoshida"

Notes

Diagonal elements are the standard realized variance of each asset (computed from their own tick data, so no synchronization needed).

K > 2 assets: implemented as O(K^2) bivariate calls. The MATLAB mfe-toolbox has a TODO here for the general case — we implement it.

Source code in src/mfe/realized/covariance.py

def realized_hayashi_yoshida(
    prices: list[FloatArray],
    times: list[FloatArray],
) -> RealizedCovarianceResult:
    """
    Hayashi-Yoshida realized covariance for K non-synchronously observed assets.

    Hayashi, T. & Yoshida, N. (2005): "On Covariance Estimation of
    Non-Synchronously Observed Diffusion Processes", Bernoulli.

    Parameters
    ----------
    prices : list of K price arrays (lengths can differ)
    times  : list of K timestamp arrays

    Returns
    -------
    RealizedCovarianceResult with (K, K) covariance matrix.
        .method = "hayashi-yoshida"

    Notes
    -----
    Diagonal elements are the standard realized variance of each asset
    (computed from their own tick data, so no synchronization needed).

    K > 2 assets: implemented as O(K^2) bivariate calls. The MATLAB
    mfe-toolbox has a TODO here for the general case — we implement it.
    """
    K = len(prices)
    if K < 2:
        raise ValueError("Need at least 2 assets.")
    if len(times) != K:
        raise ValueError(f"len(times)={len(times)} != len(prices)={K}")

    cov = np.zeros((K, K), dtype=np.float64)

    # Diagonal: realized variance from each asset's own tick data
    for k in range(K):
        r_k = np.diff(np.log(np.asarray(prices[k], dtype=np.float64)))
        cov[k, k] = float(np.sum(r_k ** 2))

    # Off-diagonal: HY estimator for each pair
    for i in range(K):
        for j in range(i + 1, K):
            hy = _hy_bivariate(
                np.asarray(prices[i], dtype=np.float64),
                np.asarray(times[i], dtype=np.float64),
                np.asarray(prices[j], dtype=np.float64),
                np.asarray(times[j], dtype=np.float64),
            )
            cov[i, j] = hy
            cov[j, i] = hy

    n_returns = min(len(p) - 1 for p in prices)

    return RealizedCovarianceResult(
        cov=cov,
        method="hayashi-yoshida",
        n_assets=K,
        n_returns=n_returns,
    )

realized_covariance_refresh_time

realized_covariance_refresh_time(prices: list[FloatArray], times: list[FloatArray]) -> RealizedCovarianceResult

Realized covariance using refresh-time synchronization.

Barndorff-Nielsen, Hansen, Lunde & Shephard (2011).

Synchronizes K asynchronous tick streams via refresh time, then computes the standard realized covariance on the synchronized returns.

Less efficient than HY (more data loss from synchronization) but gives a positive semi-definite matrix by construction.

Source code in src/mfe/realized/covariance.py

def realized_covariance_refresh_time(
    prices: list[FloatArray],
    times: list[FloatArray],
) -> RealizedCovarianceResult:
    """
    Realized covariance using refresh-time synchronization.

    Barndorff-Nielsen, Hansen, Lunde & Shephard (2011).

    Synchronizes K asynchronous tick streams via refresh time, then computes
    the standard realized covariance on the synchronized returns.

    Less efficient than HY (more data loss from synchronization) but gives a
    positive semi-definite matrix by construction.
    """
    sync_prices, sync_times = refresh_time(prices, times)
    sync_returns = np.column_stack([
        np.diff(np.log(p)) for p in sync_prices
    ])
    return realized_covariance(sync_returns)

jumps

Jump detection tests.

Barndorff-Nielsen & Shephard (2006): "Econometrics of Testing for Jumps in Financial Economics Using Bipower Variation", JFEC.

Also includes the ratio-based test and the min/med variance jump test.

bns_jump_test

bns_jump_test(returns: FloatArray, alpha: float = 0.05) -> JumpTestResult

Barndorff-Nielsen & Shephard (2006) jump test based on the ratio RV/BPV.

Z = sqrt(n) * (RV/BPV - 1) / sqrt(omega_hat)

Under the null of no jumps, Z -> N(0, 1).

Parameters:

Name	Type	Description	Default
returns	(M,) log-return array		required
alpha	float		0.05

Returns:

Type	Description
JumpTestResult

Source code in src/mfe/realized/jumps.py

def bns_jump_test(
    returns: FloatArray,
    alpha: float = 0.05,
) -> JumpTestResult:
    """
    Barndorff-Nielsen & Shephard (2006) jump test based on the ratio RV/BPV.

    Z = sqrt(n) * (RV/BPV - 1) / sqrt(omega_hat)

    Under the null of no jumps, Z -> N(0, 1).

    Parameters
    ----------
    returns : (M,) log-return array
    alpha   : significance level

    Returns
    -------
    JumpTestResult
    """
    r = np.asarray(returns, dtype=np.float64)
    n = len(r)

    rv = realized_variance(r).value
    bpv = realized_bipower_variation(r).value
    tpq = realized_tripower_quarticity(r).value

    # Consistent estimate of asymptotic variance (BN-S 2006, Theorem 1)
    # omega = (pi^2/4 + pi - 5) * max(1, TQ/BPV^2)
    pi = np.pi
    omega_hat = (pi ** 2 / 4 + pi - 5) * max(1.0, tpq / max(bpv ** 2, 1e-300))

    z_stat = float(np.sqrt(n) * (rv / max(bpv, 1e-300) - 1) / np.sqrt(max(omega_hat, 1e-300)))
    p_val = float(2 * (1 - stats.norm.cdf(abs(z_stat))))

    jump_var = max(rv - bpv, 0.0)
    return JumpTestResult(
        statistic=z_stat,
        p_value=p_val,
        jump_variation=jump_var,
        continuous_variation=bpv,
        total_variation=rv,
        significant=(p_val < alpha),
    )

kernel

Realized kernel estimator.

Barndorff-Nielsen, Hansen, Lunde & Shephard (2008): "Designing Realized Kernels to Measure the Ex-Post Variation of Equity Prices in the Presence of Noise", Econometrica.

The estimator is: RK = sum_{h=-H}^{H} k(h/(H+1)) * gamma_h where gamma_h = sum_{t>|h|} r_t * r_{t-|h|} (autocovariance of returns).

Key design decisions vs. the MATLAB version: - parameter validation is fully separated from the hot path - bandwidth selection is a standalone function (easily unit-testable) - the inner autocovariance loop is in a Cython extension (_core.pyx); if unavailable we fall back to the numpy path here

select_bandwidth

select_bandwidth(returns: FloatArray, kernel_type: KernelType = KernelType.PARZEN, noise_variance: float | None = None, iq_lower_bound: float | None = None) -> int

Optimal bandwidth H for the realized kernel.

H* = c_star * xi^{4/5} * n^{3/5}

where xi = noise_variance / sqrt(IQ), and c_star depends on the kernel.

Source code in src/mfe/realized/kernel.py

def select_bandwidth(
    returns: FloatArray,
    kernel_type: KernelType = KernelType.PARZEN,
    noise_variance: float | None = None,
    iq_lower_bound: float | None = None,
) -> int:
    """
    Optimal bandwidth H for the realized kernel.

    H* = c_star * xi^{4/5} * n^{3/5}

    where xi = noise_variance / sqrt(IQ), and c_star depends on the kernel.
    """
    r = np.asarray(returns, dtype=np.float64)
    n = len(r)

    if noise_variance is None:
        noise_variance = estimate_noise_variance(r)

    # Lower bound for IQ: use tripower variation as proxy
    if iq_lower_bound is None:
        from mfe.realized.quarticity import realized_quarticity
        iq_lower_bound = realized_quarticity(r).value

    # c_star depends on kernel (from Table 1 of BNHLS 2009)
    c_star_map = {
        KernelType.PARZEN: 3.51,
        KernelType.BARTLETT: 2.16,
        KernelType.TUKEY_HANNING: 3.68,
        KernelType.CUBIC: 3.71,
        KernelType.EPANECHNIKOV: 3.28,
        KernelType.FLAT_TOP: 2.78,
    }
    c_star = c_star_map.get(kernel_type, 3.51)

    xi = noise_variance / max(iq_lower_bound ** 0.5, 1e-30)
    H = max(1, int(np.round(c_star * (xi ** 0.4) * (n ** 0.6))))

    return H

realized_kernel

realized_kernel(returns: FloatArray, kernel_type: KernelType = KernelType.PARZEN, bandwidth: int | None = None, jitter: bool = True) -> RealizedKernelResult

Realized kernel estimator for quadratic variation.

Parameters:

Name	Type	Description	Default
returns	FloatArray		required
kernel_type	KernelType		PARZEN
bandwidth	int | None		None
jitter	bool	as in the BNHLS paper; adds a small fraction of the noise variance to handle the boundary bias	True

Returns:

Type	Description
RealizedKernelResult

Source code in src/mfe/realized/kernel.py

def realized_kernel(
    returns: FloatArray,
    kernel_type: KernelType = KernelType.PARZEN,
    bandwidth: int | None = None,
    jitter: bool = True,
) -> RealizedKernelResult:
    """
    Realized kernel estimator for quadratic variation.

    Parameters
    ----------
    returns      : (M,) log-return array (already filtered/sampled)
    kernel_type  : which kernel weight function to use
    bandwidth    : H; if None, uses automatic selector
    jitter       : if True, apply end-point jittering (noise correction)
                   as in the BNHLS paper; adds a small fraction of the
                   noise variance to handle the boundary bias

    Returns
    -------
    RealizedKernelResult
    """
    r = np.asarray(returns, dtype=np.float64)
    n = len(r)

    # Noise variance (needed for bandwidth selection and jitter)
    noise_var = estimate_noise_variance(r)

    if bandwidth is None:
        H = select_bandwidth(r, kernel_type=kernel_type, noise_variance=noise_var)
    else:
        H = int(bandwidth)

    if H >= n:
        warnings.warn(
            f"Bandwidth H={H} >= n={n}; clipping to n//2.",
            RuntimeWarning,
            stacklevel=2,
        )
        H = n // 2

    # Compute autocovariances
    if _HAS_CYTHON:
        gamma = _acov_fast(r, H)
    else:
        gamma = _autocovariance_numpy(r, H)

    # Kernel weights
    w = _kernel_weights(kernel_type, H)

    # RK = gamma_0 + 2 * sum_{h=1}^{H} k(h/(H+1)) * gamma_h
    rk = gamma[0] + 2.0 * float(w[1:] @ gamma[1:])

    # End-point (jitter) correction: adjusts for noise at boundaries
    # See BNHLS eq. (29); adds 2 * noise_var
    rk_adjusted = rk - 2.0 * noise_var if jitter else rk

    rk_adjusted = max(rk_adjusted, 0.0)  # enforce positivity

    return RealizedKernelResult(
        rk=rk,
        rk_adjusted=rk_adjusted,
        bandwidth=H,
        noise_variance=noise_var,
        iq_lower_bound=0.0,  # populated by caller if needed
        kernel_type=kernel_type,
        n_returns=n,
    )

multivariate_kernel

Multivariate Realized Kernel.

Barndorff-Nielsen, Hansen, Lunde & Shephard (2011): "Multivariate Realised Kernels: Consistent Positive Semi-Definite Estimators of the Covariation of Equity Prices with Noise and Non-Synchronous Trading", Journal of Econometrics, 162(2), 149-169.

The multivariate realized kernel generalizes the univariate realized kernel to estimate the entire (K, K) covariance matrix simultaneously from synchronised returns (after refresh-time or calendar-time sampling).

RK_{ij} = sum_{h=-H}^{H} k(h/(H+1)) * Gamma_{h,ij}

where Gamma_{h,ij} = sum_{t>|h|} r_{i,t} * r_{j,t-|h|} is the cross-autocovariance.

The key property: the full (K, K) matrix is positive semi-definite by construction because all kernels are symmetric and the weight function k satisfies k(0)=1 and the matrix {Gamma_h} has a positive-semidefinite kernel-weighted combination.

Contrast with pairwise HY: HY applied to each (i,j) pair is not guaranteed to give a PSD matrix. The multivariate realized kernel IS PSD.

Computational note: The inner loop is O(H * T * K^2). For K=10, T=50K, H=30: ~1.5B ops. The Cython _core.pyx extension covers the K=2 case analytically; for general K we use numpy's matmul broadcasting.

realized_multivariate_kernel

realized_multivariate_kernel(returns: FloatArray, kernel_type: KernelType = KernelType.PARZEN, bandwidth: int | None = None, jitter: bool = True) -> MultivariateKernelResult

Multivariate realized kernel for the full (K, K) covariance matrix.

Parameters:

Name	Type	Description	Default
returns	FloatArray		required
kernel_type	kernel weight function		PARZEN
bandwidth	int | None		None
jitter	bool	2 * diag(noise_var) from the diagonal)	True

Returns:

Type	Description
MultivariateKernelResult	.rk — raw (K, K) realized kernel matrix .rk_adjusted — noise-corrected version (PSD enforced)

Source code in src/mfe/realized/multivariate_kernel.py

def realized_multivariate_kernel(
    returns: FloatArray,
    kernel_type: KernelType = KernelType.PARZEN,
    bandwidth: int | None = None,
    jitter: bool = True,
) -> MultivariateKernelResult:
    """
    Multivariate realized kernel for the full (K, K) covariance matrix.

    Parameters
    ----------
    returns     : (T, K) synchronized log-return matrix
    kernel_type : kernel weight function
    bandwidth   : H; if None, uses average of per-asset optimal bandwidths
    jitter      : if True, apply end-point noise correction (subtract
                  2 * diag(noise_var) from the diagonal)

    Returns
    -------
    MultivariateKernelResult
        .rk          — raw (K, K) realized kernel matrix
        .rk_adjusted — noise-corrected version (PSD enforced)
    """
    R = np.asarray(returns, dtype=np.float64)
    if R.ndim == 1:
        raise ValueError("returns must be (T, K) with K >= 2 for multivariate kernel")
    T, K = R.shape

    # Per-asset noise variances and optimal bandwidths
    noise_vars = np.array([estimate_noise_variance(R[:, k]) for k in range(K)])

    if bandwidth is None:
        H_per_asset = [select_bandwidth(R[:, k], kernel_type=kernel_type,
                                        noise_variance=noise_vars[k]) for k in range(K)]
        H = max(1, int(np.round(np.mean(H_per_asset))))
    else:
        H = int(bandwidth)

    if H >= T:
        warnings.warn(f"Bandwidth H={H} >= T={T}; clipping to T//2.", RuntimeWarning, stacklevel=2)
        H = T // 2

    # Kernel weights
    w = _kernel_weights(kernel_type, H)   # (H+1,)

    # Cross-autocovariance matrices Gamma_h = R[h:].T @ R[:T-h]  for h=0..H
    # Stack: (H+1, K, K) then weight-sum
    RK = np.zeros((K, K), dtype=np.float64)

    # h = 0: R.T @ R
    RK += w[0] * (R.T @ R)

    # h = 1..H: w[h] * (Gamma_h + Gamma_h.T) — symmetric kernel
    for h in range(1, H + 1):
        Gamma_h = R[h:].T @ R[:T - h]   # (K, K)
        RK += w[h] * (Gamma_h + Gamma_h.T)

    # Jitter correction: subtract 2 * diag(noise_var) per asset
    if jitter:
        noise_correction = np.diag(2.0 * noise_vars)
        rk_adj = RK - noise_correction
    else:
        rk_adj = RK.copy()

    # Enforce PSD: project onto cone of PSD matrices
    eigvals, eigvecs = np.linalg.eigh(rk_adj)
    if np.any(eigvals < 0):
        eigvals_clipped = np.maximum(eigvals, 0.0)
        rk_adj = eigvecs @ np.diag(eigvals_clipped) @ eigvecs.T

    return MultivariateKernelResult(
        rk=RK,
        rk_adjusted=rk_adj,
        bandwidth=H,
        noise_variances=noise_vars,
        kernel_type=kernel_type,
        n_returns=T,
        n_vars=K,
    )

noise

Microstructure noise variance estimation.

Two approaches: 1. Bandi & Russell (2006): noise_var = -0.5 * mean(r_t * r_{t-1}) 2. Zhang, Mykland & Ait-Sahalia (2005): from the difference between full- frequency RV and a sub-sampled RV.

We default to the Bandi-Russell estimator as in the MATLAB mfe-toolbox.

estimate_noise_variance

estimate_noise_variance(returns: FloatArray, method: str = 'bandi-russell') -> float

Estimate the microstructure noise variance omega^2.

Parameters:

Name	Type	Description	Default
returns	(M,) log-return array at the finest available frequency		required
method	str		'bandi-russell'

Returns:

Type	Description
float — noise variance estimate (>= 0)

Source code in src/mfe/realized/noise.py

def estimate_noise_variance(
    returns: FloatArray,
    method: str = "bandi-russell",
) -> float:
    """
    Estimate the microstructure noise variance omega^2.

    Parameters
    ----------
    returns : (M,) log-return array at the finest available frequency
    method  : "bandi-russell" (default) or "zma"

    Returns
    -------
    float — noise variance estimate (>= 0)
    """
    r = np.asarray(returns, dtype=np.float64)

    if method == "bandi-russell":
        return _noise_bandi_russell(r)
    elif method == "zma":
        return _noise_zma(r)
    else:
        raise ValueError(f"Unknown noise estimation method: {method}")

quantile_var

Realized Quantile Variance.

Christensen, Oomen & Podolskij (2010): "Realised Quantile-Based Estimation of the Integrated Variance", Journal of Econometrics, 159(1), 74-98.

The realized quantile variance uses the quantile of the distribution of intra-period squared returns rather than their sum. It is jump-robust: jumps appear as extreme outliers in the distribution of r_t^2 and are downweighted by choosing an appropriate quantile below 1.

For a given quantile probability tau in (0, 1):

RQV(tau) = c(tau) * mean_{t} ( r_t^2 * 1{r_t^2 <= q_tau} ) * T / floor(tau * T)

where q_tau is the empirical tau-quantile of {r_t^2} and c(tau) is a calibration constant that ensures consistency under a pure diffusion:

c(tau) = 1 / (chi2_cdf(chi2_ppf(tau, df=1), df=1) ← same as tau for chi2(1))
       = 1 / tau (asymptotically, to leading order)

More precisely: since r_t^2 / sigma^2 ~ chi2(1) under normality, the expectation of the quantile-truncated version satisfies:

E[r_t^2 * 1{r_t^2 <= q_tau}] = sigma^2 * gamma(3/2, chi2_ppf(tau, 1)/2) / Gamma(3/2)

where gamma is the lower incomplete gamma function. We use this to calibrate.

Practical default: tau = 0.50 (median-based), which gives good jump robustness while retaining reasonable efficiency.

realized_quantile_variance

realized_quantile_variance(returns: FloatArray, tau: float = 0.5) -> RealizedQuantileVarResult

Realized quantile variance — jump-robust quadratic variation estimator.

Parameters:

Name	Type	Description	Default
returns	(T,) log-return array		required
tau	float	Lower tau → more jump-robust, less efficient. tau = 1 → equivalent to realized variance (no truncation).	0.5

Returns:

Type	Description
RealizedQuantileVarResult	.value — RQV estimate of integrated variance

Source code in src/mfe/realized/quantile_var.py

def realized_quantile_variance(
    returns: FloatArray,
    tau: float = 0.50,
) -> RealizedQuantileVarResult:
    """
    Realized quantile variance — jump-robust quadratic variation estimator.

    Parameters
    ----------
    returns : (T,) log-return array
    tau     : quantile probability in (0, 1); default 0.50 (median)
              Lower tau → more jump-robust, less efficient.
              tau = 1 → equivalent to realized variance (no truncation).

    Returns
    -------
    RealizedQuantileVarResult
        .value — RQV estimate of integrated variance
    """
    if not 0 < tau < 1:
        raise ValueError(f"tau must be in (0, 1), got {tau}")

    r = np.asarray(returns, dtype=np.float64)
    T = len(r)
    r2 = r ** 2

    # Quantile cutoff on squared returns
    q_tau = float(np.quantile(r2, tau))

    # Truncated sum: only squared returns below q_tau
    mask = r2 <= q_tau
    n_used = int(np.sum(mask))
    truncated_sum = float(np.sum(r2[mask]))

    if n_used == 0:
        return RealizedQuantileVarResult(
            value=np.nan, tau=tau, n_returns=T, n_truncated=0, quantile_cutoff=q_tau
        )

    # Raw estimate: scale so that the mean over ALL T observations is calibrated
    raw = truncated_sum / T

    # Calibration constant
    c = _calibration_constant(tau)

    rqv = raw / c if c > 0 else raw

    return RealizedQuantileVarResult(
        value=rqv,
        tau=tau,
        n_returns=T,
        n_truncated=n_used,
        quantile_cutoff=q_tau,
    )

quarticity

Realized quarticity and integrated quarticity estimators.

Used as inputs for CLT-based inference on RV and RK.

realized_quarticity

realized_quarticity(returns: FloatArray) -> RealizedResult

Realized quarticity: (n/3) * sum r_t^4

Consistent estimator of integrated quarticity IQ = int_0^1 sigma_t^4 dt.

Source code in src/mfe/realized/quarticity.py

def realized_quarticity(returns: FloatArray) -> RealizedResult:
    """
    Realized quarticity: (n/3) * sum r_t^4

    Consistent estimator of integrated quarticity IQ = int_0^1 sigma_t^4 dt.
    """
    r = np.asarray(returns, dtype=np.float64)
    n = len(r)
    rq = float((n / 3) * np.sum(r ** 4))
    return RealizedResult(value=rq, n_returns=n)

realized_tripower_quarticity

realized_tripower_quarticity(returns: FloatArray) -> RealizedResult

Tripower quarticity — robust to occasional jumps.

TPQ = n * mu_{4/3}^{-3} * mean(|r_{t-2}|^{4/3} |r_{t-1}|^{4/3} |r_t|^{4/3})

Source code in src/mfe/realized/quarticity.py

def realized_tripower_quarticity(returns: FloatArray) -> RealizedResult:
    """
    Tripower quarticity — robust to occasional jumps.

    TPQ = n * mu_{4/3}^{-3} * mean(|r_{t-2}|^{4/3} |r_{t-1}|^{4/3} |r_t|^{4/3})
    """
    r = np.asarray(returns, dtype=np.float64)
    n = len(r)
    absr = np.abs(r)
    tpq = float(
        n * (_MU43 ** -3) * np.mean(absr[:-2] ** (4 / 3) * absr[1:-1] ** (4 / 3) * absr[2:] ** (4 / 3))
    )
    return RealizedResult(value=tpq, n_returns=n)

range_

Realized range estimator.

Christensen, K. & Podolskij, M. (2007): "Realized Range-Based Estimation of Integrated Variance", Journal of Econometrics, 141(2), 323-349.

The realized range uses intra-interval high-low price ranges instead of squared returns. Under a continuous semimartingale, the range over a sub-interval [t_{j-1}, t_j] satisfies:

E[(log H_j - log L_j)^2] = 4 * log(2) * IV_j

where H_j (L_j) is the highest (lowest) price in the sub-interval. This is more efficient than squared returns due to the additional information in extreme intra-period prices.

The realized range estimator is: RR = (1 / (4 * log(2))) * sum_j (log H_j - log L_j)^2

Unlike RV it requires OHLC (open/high/low/close) data per sub-interval, which is the standard format from HFT bar data.

Also implements the normalized realized range for noise robustness.

realized_range

realized_range(high: FloatArray, low: FloatArray, open_: FloatArray | None = None, close: FloatArray | None = None) -> RealizedRangeResult

Realized range estimator from sub-interval high/low prices.

Parameters:

Name	Type	Description	Default
high	FloatArray		required
low	FloatArray		required
open_	(M,) optional — open log-price (used for Yang-Zhang correction)		None
close	(M,) optional — close log-price (used for Yang-Zhang correction)		None

Returns:

Type	Description
RealizedRangeResult	.value — realized range estimate of IV .n_intervals — M .efficiency — ~1.67 vs RV under Brownian motion

Notes

Input prices can be raw or log — the function uses log(high) - log(low) which equals log(high/low) regardless of whether inputs are already logs. If inputs are already log-prices, pass them directly; the formula is the same either way (log(e^h) - log(e^l) = h - l).

Source code in src/mfe/realized/range_.py

def realized_range(
    high: FloatArray,
    low: FloatArray,
    open_: FloatArray | None = None,
    close: FloatArray | None = None,
) -> RealizedRangeResult:
    """
    Realized range estimator from sub-interval high/low prices.

    Parameters
    ----------
    high  : (M,) highest log-price in each sub-interval (or raw price)
    low   : (M,) lowest log-price in each sub-interval
    open_ : (M,) optional — open log-price (used for Yang-Zhang correction)
    close : (M,) optional — close log-price (used for Yang-Zhang correction)

    Returns
    -------
    RealizedRangeResult
        .value — realized range estimate of IV
        .n_intervals — M
        .efficiency — ~1.67 vs RV under Brownian motion

    Notes
    -----
    Input prices can be raw or log — the function uses log(high) - log(low)
    which equals log(high/low) regardless of whether inputs are already logs.
    If inputs are already log-prices, pass them directly; the formula is
    the same either way (log(e^h) - log(e^l) = h - l).
    """
    H = np.asarray(high, dtype=np.float64)
    L = np.asarray(low, dtype=np.float64)

    if H.shape != L.shape:
        raise ValueError(f"high and low must have the same shape, got {H.shape} vs {L.shape}")

    if np.any(H < L):
        raise ValueError("high must be >= low in every interval")

    M = len(H)
    log_range_sq = (np.log(H) - np.log(L)) ** 2

    rr = float(np.sum(log_range_sq)) / _4LOG2

    # Theoretical efficiency of range vs squared-return: ~1.67 under BM
    efficiency = 1.6704  # exact: 1 / (4 log 2) * pi^2/2 * (1 - 2/pi)^{-1}; approx

    return RealizedRangeResult(value=rr, n_intervals=M, efficiency=efficiency)

realized_range_from_ticks

realized_range_from_ticks(price: FloatArray, time: FloatArray, interval_seconds: float = 300.0) -> RealizedRangeResult

Compute realized range from raw tick data by aggregating into OHLC bars.

Parameters:

Name	Type	Description	Default
price	FloatArray		required
time	FloatArray		required
interval_seconds	bar width in seconds (default 5 minutes)		300.0

Source code in src/mfe/realized/range_.py

def realized_range_from_ticks(
    price: FloatArray,
    time: FloatArray,
    interval_seconds: float = 300.0,
) -> RealizedRangeResult:
    """
    Compute realized range from raw tick data by aggregating into OHLC bars.

    Parameters
    ----------
    price            : (N,) raw tick prices
    time             : (N,) timestamps in seconds
    interval_seconds : bar width in seconds (default 5 minutes)
    """
    price = np.asarray(price, dtype=np.float64)
    time  = np.asarray(time, dtype=np.float64)

    t_start = time[0]
    t_end   = time[-1]

    # Build bar boundaries
    edges = np.arange(t_start, t_end + interval_seconds, interval_seconds)
    M = len(edges) - 1

    highs  = np.empty(M, dtype=np.float64)
    lows   = np.empty(M, dtype=np.float64)
    valid  = np.zeros(M, dtype=bool)

    for j in range(M):
        mask = (time >= edges[j]) & (time < edges[j + 1])
        if np.any(mask):
            bar_prices = price[mask]
            highs[j] = np.max(bar_prices)
            lows[j]  = np.min(bar_prices)
            valid[j] = True

    return realized_range(highs[valid], lows[valid])

sampling

Price filtering and return computation for HFT tick data.

Implements the sampling schemes from the MATLAB mfe-toolbox realized module: - Calendar-time sampling (fixed clock intervals) - Business-time sampling (fixed tick intervals) - Calendar-uniform (uniform in clock time via interpolation) - Business-uniform (uniform in tick space) - Fixed-grid sampling

All functions operate on raw tick data (price, timestamp) and return a filtered (price, time) pair ready for return computation.

price_filter

price_filter(price: FloatArray, time: FloatArray, time_type: TimeType = TimeType.SECONDS, sampling_type: SamplingType = SamplingType.CALENDAR_TIME, sampling_interval: float | int | FloatArray = 300) -> tuple[FloatArray, FloatArray]

Filter raw tick prices to a regular grid.

Parameters:

Name	Type	Description	Default
price	(N,) array of log or raw prices (function is agnostic)		required
time	FloatArray		required
time_type	how timestamps are encoded		SECONDS
sampling_type	sampling scheme		CALENDAR_TIME
sampling_interval	float | int | FloatArray	CalendarTime: seconds between samples BusinessTime: number of ticks between samples CalendarUniform / BusinessUniform: number of obs in the filtered grid Fixed: (M,) array of target times	300

Returns:

Type	Description
(filtered_price, filtered_time) — both (M,) arrays

Source code in src/mfe/realized/sampling.py

def price_filter(
    price: FloatArray,
    time: FloatArray,
    time_type: TimeType = TimeType.SECONDS,
    sampling_type: SamplingType = SamplingType.CALENDAR_TIME,
    sampling_interval: float | int | FloatArray = 300,
) -> tuple[FloatArray, FloatArray]:
    """
    Filter raw tick prices to a regular grid.

    Parameters
    ----------
    price : (N,) array of log or raw prices (function is agnostic)
    time  : (N,) timestamps in units specified by time_type
    time_type : how timestamps are encoded
    sampling_type : sampling scheme
    sampling_interval :
        - CalendarTime: seconds between samples
        - BusinessTime: number of ticks between samples
        - CalendarUniform / BusinessUniform: number of obs in the filtered grid
        - Fixed: (M,) array of target times

    Returns
    -------
    (filtered_price, filtered_time) — both (M,) arrays
    """
    price = np.asarray(price, dtype=np.float64)
    time = np.asarray(time, dtype=np.float64)

    if price.shape != time.shape:
        raise ValueError(f"price and time must have the same length, got {price.shape} vs {time.shape}")

    if sampling_type == SamplingType.CALENDAR_TIME:
        return _sample_calendar_time(price, time, float(sampling_interval))
    elif sampling_type == SamplingType.BUSINESS_TIME:
        return _sample_business_time(price, time, int(sampling_interval))
    elif sampling_type == SamplingType.CALENDAR_UNIFORM:
        return _sample_calendar_uniform(price, time, int(sampling_interval))
    elif sampling_type == SamplingType.BUSINESS_UNIFORM:
        return _sample_business_uniform(price, time, int(sampling_interval))
    elif sampling_type == SamplingType.FIXED:
        target_times = np.asarray(sampling_interval, dtype=np.float64)
        return _sample_fixed(price, time, target_times)
    else:
        raise ValueError(f"Unknown sampling_type: {sampling_type}")

returns_from_prices

returns_from_prices(price: FloatArray, log: bool = True) -> FloatArray

Compute log or simple returns from a price series.

Parameters:

Name	Type	Description	Default
price	(M,) filtered price array		required
log	bool		True

Returns:

Type	Description
(M - 1,) return array

Source code in src/mfe/realized/sampling.py

def returns_from_prices(price: FloatArray, log: bool = True) -> FloatArray:
    """
    Compute log or simple returns from a price series.

    Parameters
    ----------
    price : (M,) filtered price array
    log   : if True (default), use log-price differences

    Returns
    -------
    (M - 1,) return array
    """
    price = np.asarray(price, dtype=np.float64)
    if log:
        return np.diff(np.log(price))
    else:
        return np.diff(price) / price[:-1]

refresh_time

refresh_time(prices: list[FloatArray], times: list[FloatArray]) -> tuple[list[FloatArray], FloatArray]

Synchronize K asynchronous price series via refresh-time sampling (Barndorff-Nielsen et al. 2011).

For two assets this is O(N1 + N2) and vectorized. For K > 2 this loops over assets — TODO: Cython for K > 10.

Parameters:

Name	Type	Description	Default
prices	list of K (N_k,) price arrays		required
times	list[FloatArray]		required

Returns:

Name	Type	Description
sync_prices	list of K (M,) synchronized price arrays
sync_times	(M,) refresh times

Source code in src/mfe/realized/sampling.py

def refresh_time(
    prices: list[FloatArray],
    times: list[FloatArray],
) -> tuple[list[FloatArray], FloatArray]:
    """
    Synchronize K asynchronous price series via refresh-time sampling
    (Barndorff-Nielsen et al. 2011).

    For two assets this is O(N1 + N2) and vectorized.
    For K > 2 this loops over assets — TODO: Cython for K > 10.

    Parameters
    ----------
    prices : list of K (N_k,) price arrays
    times  : list of K (N_k,) time arrays (same units)

    Returns
    -------
    sync_prices : list of K (M,) synchronized price arrays
    sync_times  : (M,) refresh times
    """
    K = len(prices)
    if K < 2:
        raise ValueError("Need at least 2 assets for refresh-time synchronization.")

    # Initial refresh time: first time all assets have a quote
    t_start = max(t[0] for t in times)

    # Build the synchronized grid iteratively
    sync_times = []
    current_idx = [np.searchsorted(times[k], t_start, side="right") - 1 for k in range(K)]
    current_idx = [max(0, i) for i in current_idx]

    while True:
        # Current refresh time = max of current "last trade" times per asset
        t_refresh = max(times[k][current_idx[k]] for k in range(K))
        sync_times.append(t_refresh)

        # Advance each asset to the first tick at or after t_refresh
        new_idx = []
        for k in range(K):
            i = np.searchsorted(times[k], t_refresh, side="left")
            new_idx.append(i)

        # Check if we've exhausted any asset
        if any(new_idx[k] >= len(times[k]) for k in range(K)):
            break

        current_idx = new_idx

    sync_times_arr = np.array(sync_times, dtype=np.float64)
    sync_prices = []
    for k in range(K):
        idx = np.searchsorted(times[k], sync_times_arr, side="right") - 1
        idx = np.clip(idx, 0, len(prices[k]) - 1)
        sync_prices.append(prices[k][idx])

    return sync_prices, sync_times_arr

tsrv

Two-Scale Realized Variance (TSRV).

Zhang, L., Mykland, P.A. & Ait-Sahalia, Y. (2005): "A Tale of Two Time Scales: Determining Integrated Volatility With Noisy High-Frequency Data", JASA, 100(472), 1394-1411.

TSRV is a bias-corrected realized variance that is consistent under i.i.d. microstructure noise. It uses two sampling frequencies:

TSRV = (1 / (1 - n_slow/n_fast)) * (RV_slow - (n_slow/n_fast) * RV_fast)

where: RV_slow = realized variance at a slower (e.g. 5-min) frequency RV_fast = realized variance at the fastest available frequency (all ticks) n_slow = number of slow-scale returns n_fast = number of fast-scale returns

The correction removes the leading noise bias term from RV_fast.

Contrast with preaveraged RV (Jacod et al. 2009): Pre-averaging is an alternative noise-robust estimator that is also consistent and semiparametrically efficient, but TSRV is simpler and easier to interpret as a bias correction of standard RV.

Multi-scale version (MSRV) is implemented as an extension.

Reference

Zhang (2006): "Efficient Estimation of Stochastic Volatility Using Noisy Observations: A Multi-Scale Approach", Bernoulli.

tsrv

tsrv(returns: FloatArray, K: int | None = None) -> TSRVResult

Two-Scale Realized Variance.

Parameters:

Name	Type	Description	Default
returns	(N,) all-tick log-returns (finest available frequency)		required
K	int | None	If None, uses the optimal K* from Zhang et al. (2005): K* = (N/12)^{1/3} * (sigma_omega^2 / IQ)^{2/3} (approx N^{1/3})	None

Returns:

Type	Description
TSRVResult	.tsrv — the TSRV estimate of integrated variance .rv_fast — all-tick RV (noise-contaminated) .rv_slow — sub-sampled RV at scale K

Source code in src/mfe/realized/tsrv.py

def tsrv(
    returns: FloatArray,
    K: int | None = None,
) -> TSRVResult:
    """
    Two-Scale Realized Variance.

    Parameters
    ----------
    returns : (N,) all-tick log-returns (finest available frequency)
    K       : sub-sampling scale for the slow-scale estimator.
              If None, uses the optimal K* from Zhang et al. (2005):
              K* = (N/12)^{1/3} * (sigma_omega^2 / IQ)^{2/3}  (approx N^{1/3})

    Returns
    -------
    TSRVResult
        .tsrv  — the TSRV estimate of integrated variance
        .rv_fast — all-tick RV (noise-contaminated)
        .rv_slow — sub-sampled RV at scale K
    """
    r = np.asarray(returns, dtype=np.float64)
    N = len(r)

    if K is None:
        # Approximate optimal K: floor(N^{1/3}) from the asymptotic formula
        K = max(2, int(np.floor(N ** (1 / 3))))

    # Fast-scale RV: all ticks
    rv_fast = float(np.sum(r ** 2))
    n_fast = N

    # Slow-scale RV: average of K sub-grids at spacing K
    # Sub-grid s uses returns r[s], r[s+K], r[s+2K], ...
    # For each sub-grid, compute RV and average over s = 0..K-1
    rv_grids = np.empty(K, dtype=np.float64)
    n_slow_total = 0
    for s in range(K):
        r_sub = r[s::K]
        rv_grids[s] = float(np.sum(r_sub ** 2))
        n_slow_total += len(r_sub)

    rv_slow = float(np.mean(rv_grids))
    n_slow = n_slow_total // K  # average grid size

    # TSRV bias correction
    # RV_fast ≈ IQ_true + 2*N*omega^2  (noise bias)
    # RV_slow ≈ IQ_true + 2*n_slow*omega^2
    # TSRV = RV_slow - (n_slow/N) * RV_fast  (zeroes out the 2*omega^2 term)
    adj = n_slow / n_fast
    tsrv_raw = rv_slow - adj * rv_fast
    # Scale correction: 1/(1 - n_slow/N)
    scale = 1.0 / (1.0 - adj)
    tsrv_val = scale * tsrv_raw

    # Implied noise variance: omega^2 = (RV_fast - IQ_approx) / (2*N)
    # Use TSRV as the IQ approximation
    noise_var = max(0.0, (rv_fast - tsrv_val) / (2 * N))

    return TSRVResult(
        tsrv=tsrv_val,
        rv_fast=rv_fast,
        rv_slow=rv_slow,
        noise_variance=noise_var,
        n_fast=n_fast,
        n_slow=n_slow,
        K=K,
    )

msrv

msrv(returns: FloatArray, n_scales: int | None = None) -> MSRVResult

Multi-Scale Realized Variance (MSRV) — Zhang (2006).

Combines J sub-sampled RVs with optimally chosen weights to achieve the best rate of convergence under i.i.d. microstructure noise.

Parameters:

Name	Type	Description	Default
returns	FloatArray		required
n_scales	number of scales J; if None uses min(N^{1/2}, 30)		None

Source code in src/mfe/realized/tsrv.py

def msrv(
    returns: FloatArray,
    n_scales: int | None = None,
) -> MSRVResult:
    """
    Multi-Scale Realized Variance (MSRV) — Zhang (2006).

    Combines J sub-sampled RVs with optimally chosen weights to achieve
    the best rate of convergence under i.i.d. microstructure noise.

    Parameters
    ----------
    returns  : (N,) all-tick log-returns
    n_scales : number of scales J; if None uses min(N^{1/2}, 30)
    """
    r = np.asarray(returns, dtype=np.float64)
    N = len(r)

    J = n_scales if n_scales is not None else min(30, max(2, int(np.sqrt(N))))

    # Sub-sampled RVs at scales K = 1, 2, ..., J
    rv_scales = np.empty(J, dtype=np.float64)
    n_returns = np.empty(J, dtype=np.float64)

    for j, K in enumerate(range(1, J + 1)):
        rv_grids = [float(np.sum(r[s::K] ** 2)) for s in range(K)]
        rv_scales[j] = float(np.mean(rv_grids))
        n_returns[j] = N / K

    # Optimal weights from Zhang (2006) eq. (3.13)
    # w_K = 12 * K * (K - 1/2) * (1/n_K) / (J(J+1)(2J+1))
    K_arr = np.arange(1, J + 1, dtype=float)
    denom = J * (J + 1) * (2 * J + 1)
    weights = 12 * K_arr * (K_arr - 0.5) / (n_returns * denom)

    msrv_val = float(weights @ rv_scales)
    rv_fast = float(np.sum(r ** 2))

    return MSRVResult(
        msrv=msrv_val,
        rv_fast=rv_fast,
        n_scales=J,
        weights=weights,
        rv_per_scale=rv_scales,
    )

variance

Realized variance estimators.

All functions take log-returns (not prices) as input. The caller is responsible for sampling/filtering via realized.sampling.

References

Andersen & Bollerslev (1998) — realized variance Barndorff-Nielsen & Shephard (2004) — bipower variation Barndorff-Nielsen et al. (2008) — pre-averaged bipower variation Christensen & Podolskij (2007) — realized range Andersen, Dobrev & Schaumburg (2012) — realized min/med variance

realized_variance

realized_variance(returns: FloatArray, subsamples: int = 1) -> RealizedResult

Standard realized variance: sum of squared returns.

Parameters:

Name	Type	Description	Default
returns	FloatArray		required
subsamples	number of sub-grids for sub-sampling bias correction		1

Returns:

Type	Description
RealizedResult with .value = RV and .subsampled_value = sub-sampled RV

Source code in src/mfe/realized/variance.py

def realized_variance(
    returns: FloatArray,
    subsamples: int = 1,
) -> RealizedResult:
    """
    Standard realized variance: sum of squared returns.

    Parameters
    ----------
    returns    : (M,) log-return array
    subsamples : number of sub-grids for sub-sampling bias correction

    Returns
    -------
    RealizedResult with .value = RV and .subsampled_value = sub-sampled RV
    """
    r = np.asarray(returns, dtype=np.float64)
    rv = float(np.sum(r ** 2))

    rv_ss = None
    if subsamples > 1:
        ss_rvs = []
        for s in range(subsamples):
            r_sub = r[s::subsamples]
            ss_rvs.append(float(np.sum(r_sub ** 2)))
        rv_ss = float(np.mean(ss_rvs)) * subsamples  # scale back to full-sample

    return RealizedResult(
        value=rv,
        subsampled_value=rv_ss,
        n_returns=len(r),
    )

realized_bipower_variation

realized_bipower_variation(returns: FloatArray, skip: int = 0, subsamples: int = 1) -> RealizedResult

Realized bipower variation (BPV) with optional skip-k extension.

BPV = mu_1^{-2} * sum_{t=skip+2}^{T} |r_t| * |r_{t-skip-1}|

Parameters:

Name	Type	Description	Default
returns	FloatArray		required
skip	int		0
subsamples	sub-sampling replications for bias correction		1

Returns:

Type	Description
RealizedResult	.value = BPV .debiased_value = BPV * m/(m - skip - 1) where m = number of returns used

Source code in src/mfe/realized/variance.py

def realized_bipower_variation(
    returns: FloatArray,
    skip: int = 0,
    subsamples: int = 1,
) -> RealizedResult:
    """
    Realized bipower variation (BPV) with optional skip-k extension.

    BPV = mu_1^{-2} * sum_{t=skip+2}^{T} |r_t| * |r_{t-skip-1}|

    Parameters
    ----------
    returns    : (M,) log-return array
    skip       : number of returns to skip between the two absolute returns
    subsamples : sub-sampling replications for bias correction

    Returns
    -------
    RealizedResult
        .value           = BPV
        .debiased_value  = BPV * m/(m - skip - 1) where m = number of returns used
    """
    r = np.asarray(returns, dtype=np.float64)
    bpv, m = _bpv_core(r, skip)

    debiased = bpv * m / (m - skip - 1) if m > skip + 1 else np.nan

    bpv_ss = None
    if subsamples > 1:
        ss_vals = []
        for s in range(subsamples):
            r_sub = r[s::subsamples]
            val, m_sub = _bpv_core(r_sub, skip)
            ss_vals.append(val * subsamples)
        bpv_ss = float(np.mean(ss_vals))

    return RealizedResult(
        value=bpv,
        subsampled_value=bpv_ss,
        debiased_value=float(debiased),
        n_returns=len(r),
    )

realized_med_variance

realized_med_variance(returns: FloatArray) -> RealizedResult

Median realized variance: robust to jumps.

MedRV = (pi / (6 - 4*sqrt(3) + pi)) * (M/(M-2)) * sum_{t=2}^{T-1} median(|r_{t-1}|, |r_t|, |r_{t+1}|)^2

Vectorized: uses np.partition (O(N), not O(N log N)) to find the median of each triplet without sorting. ~3x faster than the column_stack approach.

Source code in src/mfe/realized/variance.py

def realized_med_variance(returns: FloatArray) -> RealizedResult:
    """
    Median realized variance: robust to jumps.

    MedRV = (pi / (6 - 4*sqrt(3) + pi)) * (M/(M-2)) *
            sum_{t=2}^{T-1} median(|r_{t-1}|, |r_t|, |r_{t+1}|)^2

    Vectorized: uses np.partition (O(N), not O(N log N)) to find the median
    of each triplet without sorting. ~3x faster than the column_stack approach.
    """
    r = np.asarray(returns, dtype=np.float64)
    M = len(r)
    absr = np.abs(r)

    a0 = absr[:-2]
    a1 = absr[1:-1]
    a2 = absr[2:]

    if _HAS_CYTHON:
        raw_sum = float(_medvar_cy(np.ascontiguousarray(absr, dtype=np.float64)))
    else:
        triplets = np.stack([a0, a1, a2], axis=1)
        partitioned = np.partition(triplets, kth=1, axis=1)
        raw_sum = float(np.sum(partitioned[:, 1] ** 2))

    pi = np.pi
    scale = (pi / (6 - 4 * np.sqrt(3) + pi)) * (M / (M - 2))
    med_rv = float(scale * raw_sum)

    return RealizedResult(value=med_rv, n_returns=M)

realized_min_variance

realized_min_variance(returns: FloatArray) -> RealizedResult

Min realized variance: minimum of adjacent pairs of squared returns.

MinRV = (pi / (pi - 2)) * (M/(M-1)) * sum_{t=1}^{T-1} min(|r_t|, |r_{t+1}|)^2

Source code in src/mfe/realized/variance.py

def realized_min_variance(returns: FloatArray) -> RealizedResult:
    """
    Min realized variance: minimum of adjacent pairs of squared returns.

    MinRV = (pi / (pi - 2)) * (M/(M-1)) * sum_{t=1}^{T-1} min(|r_t|, |r_{t+1}|)^2
    """
    r = np.asarray(returns, dtype=np.float64)
    M = len(r)
    absr = np.abs(r)

    pairs = np.column_stack([absr[:-1], absr[1:]])
    min_sq = np.min(pairs, axis=1) ** 2
    pi = np.pi
    scale = (pi / (pi - 2)) * (M / (M - 1))
    min_rv = float(scale * np.sum(min_sq))

    return RealizedResult(value=min_rv, n_returns=M)

realized_preaveraged_variance

realized_preaveraged_variance(returns: FloatArray, theta: float = 0.8) -> RealizedResult

Pre-averaged realized variance (Jacod et al. 2009).

Uses a linear pre-averaging kernel g(x) = min(x, 1-x) with block size k_n = floor(theta * sqrt(n)).

This estimator is consistent even under microstructure noise.

Parameters:

Name	Type	Description	Default
returns	(M,) log-return array		required
theta	float		0.8

Source code in src/mfe/realized/variance.py

def realized_preaveraged_variance(
    returns: FloatArray,
    theta: float = 0.8,
) -> RealizedResult:
    """
    Pre-averaged realized variance (Jacod et al. 2009).

    Uses a linear pre-averaging kernel g(x) = min(x, 1-x) with block size
    k_n = floor(theta * sqrt(n)).

    This estimator is consistent even under microstructure noise.

    Parameters
    ----------
    returns : (M,) log-return array
    theta   : tuning parameter controlling block size (default 0.8)
    """
    r = np.asarray(returns, dtype=np.float64)
    n = len(r)
    kn = max(2, int(np.floor(theta * np.sqrt(n))))

    # g(x) = min(x, 1-x) evaluated at x = i/kn for i=1..kn-1
    i_vals = np.arange(1, kn)
    g = np.minimum(i_vals / kn, 1 - i_vals / kn)
    g_sq_sum = float(np.sum(g ** 2))
    psi2 = g_sq_sum / kn  # psi_2 normalization constant

    # Pre-average
    pre_avg = np.zeros(n - kn + 1, dtype=np.float64)
    for j in range(kn - 1):
        pre_avg[: n - kn + 1] += g[j] * r[j : n - kn + 1 + j]

    pv = float(np.sum(pre_avg ** 2)) / (kn * psi2)

    # Noise bias correction (uses realized variance at fine scale)
    rv_fine = float(np.sum(r ** 2))
    bias = (kn / 2) * psi2 * rv_fine
    pv_corrected = pv - bias / kn

    return RealizedResult(
        value=pv_corrected,
        n_returns=n,
        diagnostics={"kn": kn, "theta": theta, "psi2": psi2},
    )

realized_semivariance

realized_semivariance(returns: FloatArray) -> tuple[RealizedResult, RealizedResult]

Decompose RV into positive and negative semivariance.

RS+ = sum_{r > 0} r^2, RS- = sum_{r < 0} r^2

Returns (rs_pos, rs_neg).

Source code in src/mfe/realized/variance.py

def realized_semivariance(
    returns: FloatArray,
) -> tuple[RealizedResult, RealizedResult]:
    """
    Decompose RV into positive and negative semivariance.

    RS+ = sum_{r > 0} r^2,  RS- = sum_{r < 0} r^2

    Returns (rs_pos, rs_neg).
    """
    r = np.asarray(returns, dtype=np.float64)
    rs_pos = float(np.sum(r[r > 0] ** 2))
    rs_neg = float(np.sum(r[r < 0] ** 2))
    return (
        RealizedResult(value=rs_pos, n_returns=int(np.sum(r > 0))),
        RealizedResult(value=rs_neg, n_returns=int(np.sum(r < 0))),
    )