Volatility Modeling (wraquant.vol)

The volatility module provides 28+ functions for measuring, modeling, and forecasting financial volatility – the most critical input to risk management, option pricing, and portfolio construction.

Four families of volatility tools:

  1. Realized volatility estimators – non-parametric measures from OHLC data (close-to-close, Parkinson, Garman-Klass, Rogers-Satchell, Yang-Zhang)

  2. GARCH family – parametric conditional volatility models for forecasting (GARCH, EGARCH, GJR, FIGARCH, HARCH, APARCH, DCC)

  3. Stochastic and self-exciting vol – latent process models (stochastic vol, Hawkes processes, Gaussian mixture vol)

  4. Implied vol and VRP – market-derived measures (SVI surface, variance risk premium)

Quick Example

from wraquant.vol import yang_zhang, garch_fit, egarch_fit, garch_forecast

# Best general-purpose realized vol from OHLC data
rv = yang_zhang(open, high, low, close, window=21)
print(f"Current realized vol: {rv.iloc[-1]:.4f}")

# Fit GARCH(1,1) for conditional volatility
result = garch_fit(returns)
print(f"alpha={result['alpha']:.4f}, beta={result['beta']:.4f}")
print(f"Persistence: {result['alpha'] + result['beta']:.4f}")

# EGARCH for leverage effect
egarch = egarch_fit(returns)
print(f"Leverage parameter: {egarch['gamma']:.4f}")

# 10-day ahead forecast
forecast = garch_forecast(result, horizon=10)
print(f"1-day vol: {forecast['forecasts'][0]:.4f}")

Model Selection

from wraquant.vol import garch_model_selection

# Automatic selection across GARCH variants
selection = garch_model_selection(returns)
print(f"Best model: {selection['best_model']}")
for m in selection['rankings']:
    print(f"  {m['name']:<12} BIC={m['bic']:.2f}")

Multivariate Volatility (DCC)

from wraquant.vol import dcc_fit

dcc = dcc_fit(multi_returns)
print(f"Time-varying correlation:\n{dcc['correlations'].iloc[-1]}")

See also

API Reference

Volatility modeling and forecasting.

This module provides a complete toolkit for measuring, modeling, and forecasting financial volatility – the most critical input to risk management, option pricing, and portfolio construction.

Key concepts

Volatility is the degree of variation in asset returns over time. Unlike prices, volatility is not directly observable and must be estimated. The methods here fall into four families:

  1. Realized volatility estimators – non-parametric measures computed from observed prices. Use these for ex-post volatility measurement.

    • realized_volatility – close-to-close standard deviation, the simplest estimator. Works well with daily data but ignores intraday information.

    • parkinson – uses high and low prices, ~5x more efficient than close-to-close for continuous diffusions.

    • garman_klass – uses open/high/low/close, the most efficient single-day estimator (~8x close-to-close).

    • rogers_satchell – handles drift (trending markets) unlike Parkinson; preferred when the asset has a non-zero mean return.

    • yang_zhang – combines overnight (open-to-close) and intraday (Rogers-Satchell) components; the best general-purpose estimator for daily OHLC data.

  2. GARCH family – parametric conditional volatility models. Use these when you need volatility forecasts and want to capture clustering (large moves beget large moves).

    • garch_fit – standard GARCH(1,1). The workhorse model; start here.

    • egarch_fit – Exponential GARCH. Captures the leverage effect (negative returns increase vol more than positive returns) without requiring positivity constraints.

    • gjr_garch_fit – GJR-GARCH (Glosten-Jagannathan-Runkle). Also captures leverage via an asymmetric threshold term. Simpler than EGARCH, often fits equity indices well.

    • figarch_fit – Fractionally Integrated GARCH. Models long-memory in volatility (slow hyperbolic decay of shocks). Use for assets where vol persistence is extremely high (e.g., FX, commodities).

    • harch_fit – Heterogeneous ARCH. Mixes multiple time horizons (daily, weekly, monthly) reflecting heterogeneous market participants.

    • garch_forecast – multi-step ahead volatility forecast from any fitted GARCH model.

    • dcc_fit – Dynamic Conditional Correlation GARCH for multivariate volatility and time-varying correlations across assets.

  3. Stochastic and self-exciting volatility – models where volatility itself is a latent random process or event-driven.

    • stochastic_vol_sv – Heston-style stochastic volatility estimated via particle filter or MCMC. Use when GARCH is too rigid, e.g., for option-implied dynamics.

    • hawkes_process – self-exciting point process for modeling volatility clustering through event arrivals (jumps, flash crashes).

    • gaussian_mixture_vol – regime-switching volatility using Gaussian mixtures. Use when you suspect distinct market regimes (calm vs. crisis).

  4. Implied volatility and variance risk premium – market-derived measures.

    • vol_surface_svi – fit a Stochastic Volatility Inspired (SVI) parameterisation to an implied volatility surface. Use for option pricing and vol arbitrage.

    • variance_risk_premium – difference between implied and realized variance. Positive VRP means options are “expensive” relative to realized moves.

Diagnostics

  • news_impact_curve – visualise the asymmetric response of conditional variance to return shocks (leverage effect).

  • volatility_persistence – measure how long volatility shocks take to decay (sum of GARCH alpha + beta).

How to choose

  • Daily OHLC data, one asset: start with yang_zhang for realized vol, then fit garch_fit for forecasts. If residuals show leverage, switch to egarch_fit or gjr_garch_fit.

  • Multi-asset portfolio: use dcc_fit for time-varying correlation and covariance.

  • Options / derivatives: use vol_surface_svi for the smile, and variance_risk_premium to gauge richness.

  • Regime-aware strategies: combine gaussian_mixture_vol with the regimes module.

References

  • Bollerslev (1986), “Generalized Autoregressive Conditional Heteroskedasticity”

  • Engle (2002), “Dynamic Conditional Correlation”

  • Gatheral (2004), “A parsimonious arbitrage-free implied volatility parameterization” (SVI)

  • Yang & Zhang (2000), “Drift Independent Volatility Estimation”

realized_volatility(returns, window=20, annualize=True, periods_per_year=252)[source]

Rolling realized volatility from return series.

Parameters:

  • returns (Series) – Return series.

  • window (int, default: 20) – Rolling window size.

  • annualize (bool, default: True) – Whether to annualize the volatility.

  • periods_per_year (int, default: 252) – Periods per year for annualization.

Return type:

Series

Returns:

Rolling realized volatility series.

parkinson(high, low, window=20, annualize=True, periods_per_year=252)[source]

Parkinson (1980) range-based volatility estimator.

More efficient than close-to-close for continuous processes.

Parameters:

  • high (Series) – High price series.

  • low (Series) – Low price series.

  • window (int, default: 20) – Rolling window size.

  • annualize (bool, default: True) – Whether to annualize.

  • periods_per_year (int, default: 252) – Periods per year.

Return type:

Series

Returns:

Parkinson volatility series.

garman_klass(open_, high, low, close, window=20, annualize=True, periods_per_year=252)[source]

Garman-Klass (1980) volatility estimator.

Uses open, high, low, close for higher efficiency than Parkinson or close-to-close.

Parameters:

  • open – Open price series.

  • high (Series) – High price series.

  • low (Series) – Low price series.

  • close (Series) – Close price series.

  • window (int, default: 20) – Rolling window size.

  • annualize (bool, default: True) – Whether to annualize.

  • periods_per_year (int, default: 252) – Periods per year.

  • open_ (Series)

Return type:

Series

Returns:

Garman-Klass volatility series.

rogers_satchell(open_, high, low, close, window=20, annualize=True, periods_per_year=252)[source]

Rogers-Satchell (1991) volatility estimator.

Handles drift in the price process, unlike Parkinson or Garman-Klass.

Parameters:

  • open – Open price series.

  • high (Series) – High price series.

  • low (Series) – Low price series.

  • close (Series) – Close price series.

  • window (int, default: 20) – Rolling window size.

  • annualize (bool, default: True) – Whether to annualize.

  • periods_per_year (int, default: 252) – Periods per year.

  • open_ (Series)

Return type:

Series

Returns:

Rogers-Satchell volatility series.

yang_zhang(open_, high, low, close, window=20, annualize=True, periods_per_year=252)[source]

Yang-Zhang (2000) volatility estimator.

Combines overnight and Rogers-Satchell estimators for minimum variance under drift and opening jumps.

Parameters:

  • open – Open price series.

  • high (Series) – High price series.

  • low (Series) – Low price series.

  • close (Series) – Close price series.

  • window (int, default: 20) – Rolling window size.

  • annualize (bool, default: True) – Whether to annualize.

  • periods_per_year (int, default: 252) – Periods per year.

  • open_ (Series)

Return type:

Series

Returns:

Yang-Zhang volatility series.

bipower_variation(returns, window=20, annualize=True, periods_per_year=252)[source]

Bipower variation – jump-robust realised volatility estimator.

The bipower variation (BPV) of Barndorff-Nielsen and Shephard (2004) estimates the continuous (diffusive) component of quadratic variation by using the product of adjacent absolute returns instead of squared returns. Because jumps affect only isolated observations, taking the product of consecutive absolute returns washes out the jump contribution asymptotically:

\[BPV_t = \frac{\pi}{2} \sum_{i=2}^{n} |r_{t,i}| \cdot |r_{t,i-1}|\]

where the \(\pi / 2\) factor corrects for the bias introduced by using absolute values of normal random variables (\(E[|z|] = \sqrt{2/\pi}\) for standard normal z).

Parameters:

  • returns (Series) – Return series (not prices).

  • window (int, default: 20) – Rolling window size. Each BPV estimate uses window observations.

  • annualize (bool, default: True) – Whether to annualize by sqrt(periods_per_year).

  • periods_per_year (int, default: 252) – Trading periods per year (252 for daily).

Return type:

Series

Returns:

Rolling bipower variation series (as volatility, i.e. the square root of the bipower variation divided by the window) with the same index as the input.

Example

>>> import pandas as pd, numpy as np
>>> rng = np.random.default_rng(42)
>>> returns = pd.Series(rng.normal(0, 0.01, 200))
>>> bpv = bipower_variation(returns, window=20)
>>> (bpv.dropna() > 0).all()
True

Notes

In the absence of jumps, BPV converges to the same value as standard realised variance. When jumps are present, BPV < RV and the difference RV - BPV estimates the jump component.

This is a rolling window version that computes BPV over a trailing window of window observations at each point.

Reference: Barndorff-Nielsen, O.E. and Shephard, N. (2004). “Power and Bipower Variation with Stochastic Volatility and Jumps.” Journal of Financial Econometrics, 2(1), 1–37.

See also

realized_volatility: Standard rolling realised volatility. jump_test_bns: Formal jump detection test using BPV.

jump_test_bns(returns, window=20, alpha=0.05)[source]

Barndorff-Nielsen & Shephard (2004) jump detection test.

Tests whether jumps are present in a return series by comparing realised variance (RV) to bipower variation (BPV). Under the null hypothesis of no jumps the two estimators are asymptotically equivalent. A significantly positive difference indicates the presence of jumps.

The test statistic is:

\[z = \frac{RV - BPV}{\sqrt{\theta \cdot QP}}\]

where QP is the realised quarticity (estimated via the quad-power quarticity estimator) and theta is a finite-sample correction factor:

\[\theta = \left(\frac{\pi^2}{4} + \pi - 5\right) \cdot \frac{1}{n}\]

The statistic is asymptotically standard normal under the null of no jumps.

Parameters:

  • returns (Series) – Return series (not prices).

  • window (int, default: 20) – Number of observations to use for the test. The test is applied to the most recent window observations. Use the full series length by setting window=len(returns).

  • alpha (float, default: 0.05) – Significance level for the jump detection decision. Default 0.05.

Returns:

  • rv (float) – Realised variance over the window.

  • bpv (float) – Bipower variation over the window.

  • jump_component (float) – max(RV - BPV, 0).

  • continuous_component (float) – BPV.

  • z_statistic (float) – Test statistic (standard normal under the null of no jumps).

  • p_value (float) – One-sided p-value.

  • jump_detected (bool) – True if p-value < alpha.

Return type:

dict[str, object]

Example

>>> import pandas as pd, numpy as np
>>> rng = np.random.default_rng(42)
>>> returns = pd.Series(rng.normal(0, 0.01, 100))
>>> result = jump_test_bns(returns, window=100)
>>> 'z_statistic' in result and 'p_value' in result
True

Notes

The test has low power against small or infrequent jumps in short samples. For daily data, windows of at least 60–100 observations are recommended.

Reference: Barndorff-Nielsen, O.E. and Shephard, N. (2004). “Power and Bipower Variation with Stochastic Volatility and Jumps.” Journal of Financial Econometrics, 2(1), 1–37.

See also

bipower_variation: The jump-robust volatility estimator. realized_volatility: Standard realised volatility.

two_scale_realized_variance(returns, window=20, n_slow=5, annualize=True, periods_per_year=252)[source]

Two-Scale Realised Variance (TSRV) – noise-robust RV estimator.

The TSRV of Zhang, Mykland, and Ait-Sahalia (2005) eliminates the bias caused by market microstructure noise (bid-ask bounce, discrete prices) by combining realised variances computed at two different sampling frequencies. The fast-scale estimator is contaminated by noise; the slow-scale estimator has less noise but more discretisation error. Their optimal linear combination cancels the noise bias:

\[TSRV = RV^{(slow)} - \frac{\bar{n}}{n} RV^{(all)}\]

where \(RV^{(slow)}\) is the average of sub-sampled RVs at a slower frequency, \(RV^{(all)}\) uses all ticks, and \(\bar{n}\) is the average sub-sample size.

Parameters:

  • returns (Series) – Return series. For intraday data, these should be tick-by-tick or high-frequency log returns.

  • window (int, default: 20) – Rolling window size.

  • n_slow (int, default: 5) – Sub-sampling factor for the slow scale. The fast scale uses every observation; the slow scale uses every n_slow-th. Larger values reduce noise but increase discretisation error. Default 5.

  • annualize (bool, default: True) – Whether to annualize.

  • periods_per_year (int, default: 252) – Trading periods per year.

Return type:

Series

Returns:

Rolling TSRV series (as volatility, i.e. square root).

Example

>>> import pandas as pd, numpy as np
>>> rng = np.random.default_rng(42)
>>> returns = pd.Series(rng.normal(0, 0.01, 200))
>>> tsrv = two_scale_realized_variance(returns, window=20)
>>> (tsrv.dropna() > 0).all()
True

Notes

TSRV is consistent for the integrated variance even in the presence of i.i.d. microstructure noise, unlike standard RV which diverges as the sampling frequency increases. It is particularly useful for intraday data sampled at very high frequencies (seconds or ticks).

Reference: Zhang, L., Mykland, P.A., and Ait-Sahalia, Y. (2005). “A Tale of Two Time Scales: Determining Integrated Volatility with Noisy High-Frequency Data.” Journal of the American Statistical Association, 100(472), 1394–1411.

See also

realized_volatility: Standard RV (no noise correction). realized_kernel: Alternative noise-robust estimator.

realized_kernel(returns, window=20, kernel='parzen', annualize=True, periods_per_year=252)[source]

Realised kernel estimator – noise-robust flat-top kernel RV.

The realised kernel of Barndorff-Nielsen, Hansen, Lunde, and Shephard (2008) uses a kernel-weighted sum of autocovariances of returns to produce a consistent, non-negative estimator of integrated variance in the presence of market microstructure noise:

\[RK = \sum_{h=-H}^{H} k\!\left(\frac{h}{H+1}\right) \hat{\gamma}(h)\]

where \(\hat{\gamma}(h) = \sum_i r_i r_{i-h}\) is the sample autocovariance at lag h, and \(k(\cdot)\) is a kernel function with \(k(0) = 1\), \(k(x) = 0\) for \(|x| > 1\).

The Parzen kernel is used by default as it guarantees non-negativity:

\[\begin{split}k(x) = \begin{cases} 1 - 6x^2 + 6|x|^3 & \text{if } |x| \le 1/2 \\ 2(1-|x|)^3 & \text{if } 1/2 < |x| \le 1 \\ 0 & \text{otherwise} \end{cases}\end{split}\]
Parameters:

  • returns (Series) – Return series.

  • window (int, default: 20) – Rolling window size.

  • kernel (str, default: 'parzen') – Kernel function. Currently supported: "parzen" (default). The Parzen kernel guarantees non-negativity.

  • annualize (bool, default: True) – Whether to annualize.

  • periods_per_year (int, default: 252) – Trading periods per year.

Return type:

Series

Returns:

Rolling realised kernel volatility series (square root of the kernel-weighted variance).

Example

>>> import pandas as pd, numpy as np
>>> rng = np.random.default_rng(42)
>>> returns = pd.Series(rng.normal(0, 0.01, 200))
>>> rk = realized_kernel(returns, window=20)
>>> (rk.dropna() > 0).all()
True

Notes

The realised kernel is rate-optimal under i.i.d. noise and achieves a convergence rate of \(n^{-1/5}\) (vs \(n^{-1/4}\) for standard RV with optimal sparse sampling).

The bandwidth H is chosen automatically as \(H = \lceil c \cdot n^{3/5} \rceil\) where c is a constant typically set to 1.

Reference: Barndorff-Nielsen, O.E., Hansen, P.R., Lunde, A., and Shephard, N. (2008). “Designing Realized Kernels to Measure the Ex-Post Variation of Equity Prices in the Presence of Noise.” Econometrica, 76(6), 1481–1536.

See also

two_scale_realized_variance: Alternative noise-robust estimator. realized_volatility: Standard RV without noise correction. bipower_variation: Jump-robust estimator.

ewma_volatility(returns, span=30, annualize=True, periods_per_year=252)[source]

Exponentially weighted moving average volatility.

The EWMA model of J.P. Morgan’s RiskMetrics (1996) assigns exponentially decaying weights to past squared returns. The decay factor lambda = 1 - 2/(span+1).

\[\sigma_t^2 = \lambda \sigma_{t-1}^2 + (1 - \lambda) r_{t-1}^2\]
Parameters:

  • returns (Series) – Return series (not prices).

  • span (int, default: 30) – EWMA span parameter. The decay factor is computed as lambda = 1 - 2/(span+1). RiskMetrics uses span ~= 73 (lambda = 0.94 for daily data). Default 30.

  • annualize (bool, default: True) – Whether to annualize the volatility by multiplying by sqrt(periods_per_year).

  • periods_per_year (int, default: 252) – Trading periods per year. Use 252 for daily, 52 for weekly, 12 for monthly.

Return type:

Series

Returns:

EWMA volatility series with same index as input.

Example

>>> import pandas as pd, numpy as np
>>> rng = np.random.default_rng(42)
>>> returns = pd.Series(rng.normal(0, 0.01, 500))
>>> vol = ewma_volatility(returns, span=30)
>>> vol.iloc[-1] > 0
True

Notes

RiskMetrics (1996) recommended lambda = 0.94 for daily data and lambda = 0.97 for monthly data. The EWMA model is a restricted IGARCH(1,1) with zero intercept.

See also

garch_fit: More flexible parametric volatility model. realized_volatility: Non-parametric rolling window estimator.

garch_fit(returns, p=1, q=1, mean='Constant', dist='normal', **kwargs)[source]

Fit a GARCH(p, q) model to a return series.

The Generalized Autoregressive Conditional Heteroskedasticity model of Bollerslev (1986) captures volatility clustering – the tendency for large moves to follow large moves. The conditional variance is:

\[\sigma_t^2 = \omega + \sum_{i=1}^{p} \alpha_i \epsilon_{t-i}^2 + \sum_{j=1}^{q} \beta_j \sigma_{t-j}^2\]
Parameters:

  • returns (Series | ndarray) – Return series (not prices). Should be percentage returns or log returns, typically in the range [-10, 10].

  • p (int, default: 1) – Order of the GARCH (lagged variance) terms. Higher values capture longer volatility memory. Default 1 is standard.

  • q (int, default: 1) – Order of the ARCH (lagged squared residual) terms. Default 1.

  • mean (str, default: 'Constant') –

    Mean model specification. Options:

    • "Constant": constant mean (most common)

    • "Zero": zero mean (for demeaned returns)

    • "ARX": autoregressive with exogenous variables

  • dist (str, default: 'normal') –

    Error distribution. Options:

    • "normal": Gaussian (fastest, least flexible)

    • "t": Student’s t (captures fat tails)

    • "skewt": Hansen’s skewed t (tails + asymmetry)

    • "ged": Generalized Error Distribution

  • **kwargs (Any) – Additional keyword arguments passed to arch.arch_model() (e.g., lags for ARX mean).

Returns:

  • model_name (str) – Model identifier, e.g. "GARCH(1,1)".

  • params (dict) – Fitted parameters (omega, alpha, beta, plus distribution params if non-normal).

  • conditional_volatility (pd.Series) – Time series of estimated conditional standard deviation. Multiply by sqrt(252) to annualize if returns are daily.

  • standardized_residuals (np.ndarray) – Residuals divided by conditional std dev. Should be ~i.i.d. if model is correct.

  • aic (float) – Akaike Information Criterion (lower = better).

  • bic (float) – Bayesian Information Criterion.

  • log_likelihood (float) – Maximized log-likelihood.

  • persistence (float) – Sum of alpha + beta. Values near 1 indicate high volatility persistence (IGARCH if exactly 1).

  • half_life (float) – Number of periods for a volatility shock to decay to half its initial impact.

  • unconditional_variance (float) – Long-run variance omega / (1 - alpha - beta).

  • ljung_box (dict) – Ljung-Box test on squared standardized residuals. If p-value > 0.05, model captures vol clustering.

  • model – The fitted arch model result object for further analysis (forecasting, simulation, etc.).

Return type:

dict[str, Any]

Example

>>> import numpy as np
>>> from wraquant.vol.models import garch_fit
>>> rng = np.random.default_rng(42)
>>> returns = rng.normal(0, 1, 1000)
>>> result = garch_fit(returns, p=1, q=1, dist="normal")
>>> 0 < result['persistence'] < 2
True

Notes

The GARCH(1,1) model is the workhorse of volatility modeling. Persistence (alpha + beta) near 1 is common for financial returns. If persistence > 1, consider IGARCH or FIGARCH.

Returns are internally scaled by 100 before fitting to improve numerical stability (the arch library convention). All outputs are rescaled back to the original units.

Reference: Bollerslev, T. (1986). “Generalized Autoregressive Conditional Heteroskedasticity.” Journal of Econometrics, 31, 307–327.

See also

egarch_fit: Asymmetric GARCH capturing leverage effect. gjr_garch_fit: Alternative asymmetric specification. garch_forecast: Multi-step ahead volatility forecasting. news_impact_curve: Visualize asymmetric shock response.

egarch_fit(returns, p=1, q=1, mean='Constant', dist='normal', **kwargs)[source]

Fit an EGARCH(p, q) model (exponential GARCH).

The EGARCH model of Nelson (1991) captures the leverage effect – negative shocks increase volatility more than positive shocks of equal magnitude. The model specifies the log of conditional variance, so positivity of variance is guaranteed without parameter constraints:

\[\ln \sigma_t^2 = \omega + \sum_{i=1}^{q} \alpha_i \left( \left|z_{t-i}\right| - E\left|z_{t-i}\right| \right) + \sum_{i=1}^{q} \gamma_i z_{t-i} + \sum_{j=1}^{p} \beta_j \ln \sigma_{t-j}^2\]

where \(z_t = \epsilon_t / \sigma_t\) are standardized residuals and \(\gamma_i\) captures the asymmetric (leverage) response.

Parameters:

  • returns (Series | ndarray) – Return series (not prices).

  • p (int, default: 1) – Order of GARCH (lagged log-variance) terms.

  • q (int, default: 1) – Order of ARCH (asymmetric innovation) terms.

  • mean (str, default: 'Constant') – Mean model specification ("Constant", "Zero", "ARX").

  • dist (str, default: 'normal') – Error distribution ("normal", "t", "skewt", "ged").

  • **kwargs (Any) – Additional keyword arguments passed to arch.arch_model().

Return type:

dict[str, Any]

Returns:

Dictionary with the same structure as garch_fit(). The params dict includes gamma terms capturing the leverage effect. A negative gamma means negative shocks increase volatility more than positive shocks.

Example

>>> import numpy as np
>>> from wraquant.vol.models import egarch_fit
>>> rng = np.random.default_rng(42)
>>> returns = rng.normal(0, 1, 1000)
>>> result = egarch_fit(returns, p=1, q=1)
>>> result['conditional_volatility'].iloc[-1] > 0
True

Notes

The EGARCH model is preferred when:

  1. Leverage effects are important (equity markets).

  2. Parameter positivity constraints are problematic.

The log specification means persistence is measured by the sum of beta coefficients. Unconditional variance may not have a closed form depending on the distribution.

Reference: Nelson, D.B. (1991). “Conditional Heteroskedasticity in Asset Returns: A New Approach.” Econometrica, 59(2), 347–370.

See also

garch_fit: Symmetric GARCH model. gjr_garch_fit: Alternative asymmetric GARCH. news_impact_curve: Compare asymmetric responses across models.

gjr_garch_fit(returns, p=1, q=1, o=1, mean='Constant', dist='normal', **kwargs)[source]

Fit a GJR-GARCH(p, o, q) model (Glosten-Jagannathan-Runkle).

The GJR-GARCH model of Glosten, Jagannathan, and Runkle (1993) augments the standard GARCH with an indicator function that activates on negative shocks, capturing the leverage effect:

\[\sigma_t^2 = \omega + \sum_{i=1}^{q} \alpha_i \epsilon_{t-i}^2 + \sum_{i=1}^{o} \gamma_i \epsilon_{t-i}^2 I(\epsilon_{t-i} < 0) + \sum_{j=1}^{p} \beta_j \sigma_{t-j}^2\]

where \(I(\cdot)\) is the indicator function. The gamma coefficient captures the extra volatility increase from negative shocks (bad news).

Parameters:

  • returns (Series | ndarray) – Return series (not prices).

  • p (int, default: 1) – Order of GARCH (lagged variance) terms.

  • q (int, default: 1) – Order of ARCH (lagged squared residual) terms.

  • o (int, default: 1) – Order of asymmetric (threshold) terms. Default 1.

  • mean (str, default: 'Constant') – Mean model specification.

  • dist (str, default: 'normal') – Error distribution.

  • **kwargs (Any) – Additional keyword arguments passed to arch.arch_model().

Return type:

dict[str, Any]

Returns:

Dictionary with the same structure as garch_fit(). The params dict includes gamma[1] capturing the asymmetric response. Positive gamma means negative shocks increase volatility more.

The persistence for GJR-GARCH is alpha + beta + 0.5 * gamma (assuming symmetric distribution).

Example

>>> import numpy as np
>>> from wraquant.vol.models import gjr_garch_fit
>>> rng = np.random.default_rng(42)
>>> returns = rng.normal(0, 1, 1000)
>>> result = gjr_garch_fit(returns, p=1, q=1)
>>> 'gamma[1]' in result['params'] or 'gamma' in str(result['params'])
True

Notes

The GJR-GARCH is equivalent to TARCH (Threshold ARCH) when applied to the variance (not standard deviation). It is the default asymmetric model in many risk management applications.

For stationarity, the persistence alpha + beta + 0.5 * gamma < 1 is required.

Reference: Glosten, L.R., Jagannathan, R., and Runkle, D.E. (1993). “On the Relation between the Expected Value and the Volatility of the Nominal Excess Return on Stocks.” Journal of Finance, 48(5), 1779–1801.

See also

egarch_fit: Log-variance based asymmetric model. garch_fit: Symmetric baseline GARCH. news_impact_curve: Compare asymmetric responses.

figarch_fit(returns, p=1, q=1, mean='Constant', dist='normal', **kwargs)[source]

Fit a FIGARCH(p, d, q) model (Fractionally Integrated GARCH).

The FIGARCH model of Baillie, Bollerslev, and Mikkelsen (1996) allows for long memory in the conditional variance process. The fractional differencing parameter d is estimated from the data and captures the slow hyperbolic decay of volatility autocorrelations that is frequently observed in financial time series.

\[\sigma_t^2 = \omega + \beta \sigma_{t-1}^2 + [1 - \beta L - (1-L)^d (1 - \alpha L)] \epsilon_t^2\]

where d in (0, 1) is the fractional integration parameter, and L is the lag operator.

Parameters:

  • returns (Series | ndarray) – Return series (not prices).

  • p (int, default: 1) – Order of GARCH terms (typically 1).

  • q (int, default: 1) – Order of ARCH terms (typically 1).

  • mean (str, default: 'Constant') – Mean model specification.

  • dist (str, default: 'normal') – Error distribution.

  • **kwargs (Any) – Additional keyword arguments passed to arch.arch_model().

Return type:

dict[str, Any]

Returns:

Dictionary with the same structure as garch_fit(). The params dict includes the fractional integration parameter d. Values of d near 0 imply short memory (like GARCH), values near 0.5 imply long memory.

Example

>>> import numpy as np
>>> from wraquant.vol.models import figarch_fit
>>> rng = np.random.default_rng(42)
>>> returns = rng.normal(0, 1, 1000)
>>> result = figarch_fit(returns, p=1, q=1)
>>> result['aic'] < 0 or result['aic'] > 0  # AIC is finite
True

Notes

FIGARCH nests GARCH (d=0) and IGARCH (d=1). Empirically, most financial return series have d in (0.3, 0.5), suggesting long memory in volatility.

The arch library implements FIGARCH via the "FIGARCH" volatility specification. The power parameter defaults to 2.0 (standard FIGARCH).

Reference: Baillie, R.T., Bollerslev, T., and Mikkelsen, H.O. (1996). “Fractionally Integrated Generalized Autoregressive Conditional Heteroskedasticity.” Journal of Econometrics, 74(1), 3–30.

See also

garch_fit: Standard GARCH without long memory. harch_fit: Heterogeneous ARCH as alternative long-memory model. volatility_persistence: Measure persistence and half-life.

aparch_fit(returns, p=1, o=1, q=1, mean='Constant', dist='normal', **kwargs)[source]

Fit an APARCH(p, o, q) model (Asymmetric Power ARCH).

The APARCH model of Ding, Granger, and Engle (1993) generalises the standard GARCH by raising the conditional standard deviation to a freely estimated power delta and introducing a leverage (rotation) parameter gamma for each asymmetric term:

\[\sigma_t^\delta = \omega + \sum_{i=1}^{q} \alpha_i \left( |\epsilon_{t-i}| - \gamma_i \epsilon_{t-i} \right)^\delta + \sum_{j=1}^{p} \beta_j \sigma_{t-j}^\delta\]

The power parameter delta is estimated from the data rather than being fixed at 2 (GARCH) or 1 (TARCH). This gives the model extra flexibility to capture the shape of the conditional distribution of volatility. When delta = 2 the model collapses to GJR-GARCH; when delta = 1 it collapses to TARCH.

Parameters:

  • returns (Series | ndarray) – Return series (not prices). Should be percentage returns or log returns.

  • p (int, default: 1) – Order of GARCH (lagged powered standard deviation) terms.

  • o (int, default: 1) – Order of asymmetric (leverage) terms. Default 1.

  • q (int, default: 1) – Order of ARCH (lagged powered absolute residual) terms.

  • mean (str, default: 'Constant') –

    Mean model specification. Options:

    • "Constant": constant mean (most common)

    • "Zero": zero mean (for demeaned returns)

    • "ARX": autoregressive with exogenous variables

  • dist (str, default: 'normal') –

    Error distribution. Options:

    • "normal": Gaussian

    • "t": Student’s t (fat tails)

    • "skewt": Hansen’s skewed t

    • "ged": Generalized Error Distribution

  • **kwargs (Any) – Additional keyword arguments passed to arch.arch_model().

Returns:

  • model_name (str) – "APARCH(p,o,q)".

  • params (dict) – Fitted parameters including omega, alpha[i], gamma[i] (leverage), beta[j], and the estimated power parameter delta (often labelled delta or embedded in the volatility process params).

  • conditional_volatility (pd.Series) – Conditional standard deviation series.

  • standardized_residuals (np.ndarray) – Model residuals divided by conditional std dev.

  • aic / bic / log_likelihood (float) – Fit quality metrics.

  • persistence (float) – Volatility persistence.

  • half_life (float) – Half-life of volatility shocks.

  • unconditional_variance (float) – Long-run variance.

  • ljung_box (dict) – Ljung-Box test on squared standardized residuals.

  • model – Fitted arch result object.

Return type:

dict[str, Any]

Example

>>> import numpy as np
>>> from wraquant.vol.models import aparch_fit
>>> rng = np.random.default_rng(42)
>>> returns = rng.normal(0, 1, 1000)
>>> result = aparch_fit(returns, p=1, o=1, q=1)
>>> result['conditional_volatility'].iloc[-1] > 0
True

Notes

The key advantage of APARCH is the estimated power parameter delta. Empirical studies (Ding et al. 1993, Laurent 2004) find that delta is often close to 1 for daily equity returns rather than 2, implying the conditional standard deviation (not variance) follows a more parsimonious dynamic.

APARCH nests at least seven other ARCH-class models as special cases, making it useful for model selection via the estimated delta.

Reference: Ding, Z., Granger, C.W.J., and Engle, R.F. (1993). “A Long Memory Property of Stock Market Returns and a New Model.” Journal of Empirical Finance, 1, 83–106.

See also

garch_fit: Standard GARCH (delta fixed at 2, no leverage). gjr_garch_fit: GJR-GARCH (delta fixed at 2, with leverage). egarch_fit: Log-variance based asymmetric model.

harch_fit(returns, lags=None, mean='Constant', dist='normal', **kwargs)[source]

Fit a HARCH model (Heterogeneous ARCH).

The HARCH model of Mueller et al. (1997) aggregates squared returns over multiple time horizons to capture the heterogeneous nature of market participants operating at different frequencies (intraday traders, daily, weekly, monthly):

\[\sigma_t^2 = \omega + \sum_{j=1}^{J} \alpha_j \left( \sum_{i=1}^{l_j} \epsilon_{t-i} \right)^2 / l_j\]

where \(l_1 < l_2 < \ldots < l_J\) are the lag lengths corresponding to different time horizons.

Parameters:

  • returns (Series | ndarray) – Return series (not prices).

  • lags (list[int] | None, default: None) – List of lag lengths representing different time horizons. Default is [1, 5, 22] corresponding to daily, weekly, and monthly horizons for daily data.

  • mean (str, default: 'Constant') – Mean model specification.

  • dist (str, default: 'normal') – Error distribution.

  • **kwargs (Any) – Additional keyword arguments passed to arch.arch_model().

Return type:

dict[str, Any]

Returns:

Dictionary with the same structure as garch_fit(). The params dict includes coefficients for each aggregation horizon.

Example

>>> import numpy as np
>>> from wraquant.vol.models import harch_fit
>>> rng = np.random.default_rng(42)
>>> returns = rng.normal(0, 1, 1000)
>>> result = harch_fit(returns, lags=[1, 5, 22])
>>> result['conditional_volatility'].iloc[-1] > 0
True

Notes

HARCH is particularly useful for modeling high-frequency data where different market participants operate on vastly different time scales. It provides a parsimonious alternative to high-order ARCH or GARCH models.

Reference: Mueller, U.A., Dacorogna, M.M., Dave, R.D., Olsen, R.B., Pictet, O.V., and von Weizsacker, J.E. (1997). “Volatilities of Different Time Resolutions.” Journal of Empirical Finance, 4(2-3), 213–239.

See also

garch_fit: Standard GARCH model. figarch_fit: Long-memory GARCH via fractional integration.

garch_forecast(returns, p=1, q=1, mean='Constant', dist='normal', horizon=10, method='analytic', simulations=1000, **kwargs)[source]

Multi-step ahead volatility forecast from a GARCH model.

Fits a GARCH(p,q) model and produces multi-step forecasts of conditional variance, with optional confidence intervals via simulation.

Parameters:

  • returns (Series | ndarray) – Return series (not prices).

  • p (int, default: 1) – GARCH lag order.

  • q (int, default: 1) – ARCH lag order.

  • mean (str, default: 'Constant') – Mean model specification.

  • dist (str, default: 'normal') – Error distribution.

  • horizon (int, default: 10) – Number of steps ahead to forecast. Default 10.

  • method (str, default: 'analytic') –

    Forecasting method. Options:

    • "analytic": Closed-form multi-step variance forecast. Fastest and exact for GARCH. Does not provide intervals.

    • "simulation": Monte Carlo simulation of future paths. Provides confidence intervals but slower.

  • simulations (int, default: 1000) – Number of simulation paths (only used when method="simulation"). Default 1000.

  • **kwargs (Any) – Additional keyword arguments passed to arch.arch_model().

Returns:

  • model_name (str) – "GARCH(p,q) forecast".

  • forecast_variance (np.ndarray) – Forecasted conditional variance for each horizon step, shape (horizon,).

  • forecast_volatility (np.ndarray) – Square root of forecast variance, shape (horizon,).

  • confidence_intervals (dict | None) – If method is "simulation", contains "lower_5" and "upper_95" arrays. None for analytic method.

  • fit_result (dict) – Full GARCH fit result from garch_fit().

Return type:

dict[str, Any]

Example

>>> import numpy as np
>>> from wraquant.vol.models import garch_forecast
>>> rng = np.random.default_rng(42)
>>> returns = rng.normal(0, 1, 1000)
>>> fc = garch_forecast(returns, horizon=5)
>>> len(fc['forecast_volatility']) == 5
True

Notes

Analytic forecasts assume the model is correctly specified and use the recursive formula for multi-step variance:

\[E[\sigma_{T+h}^2 | \mathcal{F}_T] = \omega \frac{1 - (\alpha + \beta)^h}{1 - \alpha - \beta} + (\alpha + \beta)^h \sigma_{T+1}^2\]

Simulation forecasts draw from the fitted error distribution and roll the GARCH recursion forward, providing a distribution of future paths.

See also

garch_fit: Fit the underlying GARCH model. volatility_persistence: Analyze shock decay properties.

garch_rolling_forecast(returns, window=None, horizon=1, vol='GARCH', p=1, q=1, dist='normal', refit_every=1, **kwargs)[source]

Rolling or expanding window 1-step-ahead GARCH volatility forecast.

At each step the model is fit on past data (fixed window or expanding) and a horizon-step-ahead forecast is produced. The resulting series of out-of-sample forecasts can be compared against realised volatility for back-testing purposes.

\[\hat{\sigma}_{t+1|t}^2 = \omega + \alpha \epsilon_t^2 + \beta \sigma_t^2\]

The forecast is produced at each date t using only data up to t, providing a genuinely out-of-sample volatility estimate.

Parameters:

  • returns (Series | ndarray) – Return series. If a pd.Series with a DatetimeIndex, dates are preserved in the output.

  • window (int | None, default: None) – Rolling window size (number of observations). If None an expanding window is used, starting from the first 100 observations. A fixed window is common in practice (e.g., 500 or 1000 days).

  • horizon (int, default: 1) – Forecast horizon (steps ahead). Default 1 for 1-step.

  • vol (str, default: 'GARCH') – Volatility model type. Any valid arch vol spec: "GARCH", "EGARCH", "GJR-GARCH", "APARCH", etc. Default "GARCH".

  • p (int, default: 1) – GARCH order.

  • q (int, default: 1) – ARCH order.

  • dist (str, default: 'normal') – Error distribution ("normal", "t", "skewt").

  • refit_every (int, default: 1) – Refit the model every n steps. Setting this > 1 greatly speeds up the rolling forecast by re-using the previous fit for intermediate steps. Default 1 (refit every step).

  • **kwargs (Any) – Additional keyword arguments passed to arch.arch_model().

Returns:

  • forecasts (pd.Series) – Series of forecasted conditional volatility (standard deviation), indexed by the date at which the forecast was made.

  • actuals (pd.Series) – Absolute returns as a proxy for realised volatility, aligned with forecasts.

  • dates (pd.Index | np.ndarray) – Forecast dates or integer positions.

Return type:

dict[str, Any]

Example

>>> import numpy as np, pandas as pd
>>> from wraquant.vol.models import garch_rolling_forecast
>>> rng = np.random.default_rng(42)
>>> ret = pd.Series(rng.normal(0, 0.01, 300))
>>> result = garch_rolling_forecast(ret, window=200, refit_every=10)
>>> len(result['forecasts']) > 0
True
>>> (result['forecasts'] > 0).all()
True

Notes

Setting refit_every > 1 avoids refitting the model at every single step. Between refits the most recent parameter estimates are used and the conditional variance is simply rolled forward using the GARCH recursion. This is standard practice in the literature (Hansen & Lunde 2005) and can reduce computation by an order of magnitude.

When vol="EGARCH" the arch library fits the model on log variance, so internal rescaling differs. The output is always in the same units as the input returns.

Reference: Hansen, P.R. and Lunde, A. (2005). “A Forecast Comparison of Volatility Models: Does Anything Beat a GARCH(1,1)?” Journal of Applied Econometrics, 20, 873–889.

See also

garch_forecast: Single multi-step forecast from a fitted model. garch_model_selection: Compare GARCH variants systematically.

garch_model_selection(returns, p=1, q=1, distributions=None, mean='Constant', **kwargs)[source]

Compare GARCH, EGARCH, and GJR-GARCH across error distributions.

Fits every combination of volatility specification and error distribution, then ranks the models by AIC (and reports BIC, log-likelihood, and persistence for each).

This is the standard model-selection exercise recommended by Hansen & Lunde (2005) to determine which GARCH variant and error distribution best describes a given return series.

Parameters:

  • returns (Series | ndarray) – Return series (not prices).

  • p (int, default: 1) – GARCH order (applied to all models).

  • q (int, default: 1) – ARCH order (applied to all models).

  • distributions (list[str] | None, default: None) – List of error distributions to try. Defaults to ["normal", "t", "skewt"].

  • mean (str, default: 'Constant') – Mean model specification applied to all fits.

  • **kwargs (Any) – Additional keyword arguments passed to arch.arch_model().

Returns:

  • model (str) – Model name (e.g. "EGARCH-t").

  • vol_model (str) – Volatility specification.

  • distribution (str) – Error distribution.

  • aic (float) – Akaike Information Criterion.

  • bic (float) – Bayesian Information Criterion.

  • log_likelihood (float) – Maximised log-likelihood.

  • persistence (float) – Volatility persistence.

  • num_params (int) – Number of estimated parameters.

Return type:

DataFrame

Example

>>> import numpy as np
>>> from wraquant.vol.models import garch_model_selection
>>> rng = np.random.default_rng(42)
>>> returns = rng.normal(0, 1, 1000)
>>> df = garch_model_selection(returns)
>>> 'aic' in df.columns and 'bic' in df.columns
True
>>> len(df) >= 6  # 3 models x 2+ distributions
True

Notes

Models that fail to converge are silently excluded from the results table. AIC penalises complexity less than BIC, so the best AIC and best BIC model may differ – both are reported.

Reference: Hansen, P.R. and Lunde, A. (2005). “A Forecast Comparison of Volatility Models: Does Anything Beat a GARCH(1,1)?” Journal of Applied Econometrics, 20, 873–889.

See also

garch_fit: Fit a specific GARCH model. egarch_fit: Fit EGARCH. gjr_garch_fit: Fit GJR-GARCH.

dcc_fit(returns, p=1, q=1, dist='normal')[source]

Fit a DCC-GARCH model for multivariate dynamic correlations.

The Dynamic Conditional Correlation model of Engle (2002) estimates time-varying correlations between multiple asset returns using a two-step procedure:

  1. Fit univariate GARCH(p,q) to each series.

  2. Estimate DCC parameters (a, b) via MLE on standardized residuals.

The DCC correlation dynamics are:

\[ \begin{align}\begin{aligned}Q_t = (1 - a - b) \bar{Q} + a \, z_{t-1} z_{t-1}^\top + b \, Q_{t-1}\\R_t = \text{diag}(Q_t)^{-1/2} \, Q_t \, \text{diag}(Q_t)^{-1/2}\end{aligned}\end{align} \]

where \(\bar{Q}\) is the unconditional correlation matrix of standardized residuals.

Parameters:

  • returns (DataFrame | ndarray) – Return data, either a DataFrame with k columns (one per asset) or an array of shape (T, k).

  • p (int, default: 1) – GARCH lag order for each univariate model.

  • q (int, default: 1) – ARCH lag order for each univariate model.

  • dist (str, default: 'normal') – Error distribution for univariate GARCH fits.

Returns:

  • dcc_params (dict) – DCC parameters {"a": ..., "b": ...}.

  • univariate_results (list[dict]) – Per-asset GARCH fit results from garch_fit().

  • conditional_correlations (np.ndarray) – Array of shape (T, k, k) with time-varying correlation matrices.

  • conditional_covariances (np.ndarray) – Array of shape (T, k, k) with time-varying covariance matrices.

  • conditional_volatilities (np.ndarray) – Array of shape (T, k) with per-asset conditional volatilities.

  • standardized_residuals (np.ndarray) – Array of shape (T, k).

  • qbar (np.ndarray) – Unconditional correlation matrix.

Return type:

dict[str, Any]

Example

>>> import numpy as np, pandas as pd
>>> from wraquant.vol.models import dcc_fit
>>> rng = np.random.default_rng(42)
>>> ret = pd.DataFrame({
...     'A': rng.normal(0, 0.01, 500),
...     'B': rng.normal(0, 0.015, 500),
... })
>>> result = dcc_fit(ret)
>>> result['conditional_correlations'].shape == (500, 2, 2)
True

Notes

This implementation uses the arch library for the univariate GARCH step, then applies the pure scipy/numpy DCC estimation from wraquant.risk.dcc.

The DCC model is widely used in portfolio management, value-at-risk computation, and hedging applications where time-varying correlations are important.

Reference: Engle, R. (2002). “Dynamic Conditional Correlation: A Simple Class of Multivariate Generalized Autoregressive Conditional Heteroskedasticity Models.” Journal of Business & Economic Statistics, 20(3), 339–350.

See also

garch_fit: Univariate GARCH estimation. garch_forecast: Volatility forecasting.

realized_garch(returns, realized_vol, p=1, q=1, mean='Constant', dist='normal', **kwargs)[source]

Fit a Realized GARCH model augmented with realized volatility.

The Realized GARCH model of Hansen, Huang, and Shek (2012) augments the standard GARCH with a measurement equation linking the conditional variance to a realized volatility measure, improving forecasting performance compared to standard GARCH.

This implementation uses the arch library’s ARX mean model with lagged realized volatility as an exogenous regressor, capturing the predictive information in the realized measure.

Parameters:

  • returns (Series | ndarray) – Return series (not prices).

  • realized_vol (Series | ndarray) – Corresponding realized volatility measure (must have the same length as returns). Can be any realized measure: realized variance, Parkinson, Garman-Klass, etc.

  • p (int, default: 1) – GARCH lag order.

  • q (int, default: 1) – ARCH lag order.

  • mean (str, default: 'Constant') – Mean model specification. Default "Constant" uses a constant plus the realized vol regressor.

  • dist (str, default: 'normal') – Error distribution.

  • **kwargs (Any) – Additional keyword arguments passed to arch.arch_model().

Return type:

dict[str, Any]

Returns:

Dictionary with the same structure as garch_fit(), plus:

  • realized_vol_used (np.ndarray) – The realized volatility series used in the model.

Example

>>> import numpy as np
>>> from wraquant.vol.models import realized_garch
>>> rng = np.random.default_rng(42)
>>> returns = rng.normal(0, 0.01, 500)
>>> rv = np.abs(returns).rolling(20).mean() if hasattr(returns, 'rolling') else np.convolve(np.abs(returns), np.ones(20)/20, mode='same')
>>> # result = realized_garch(returns, rv)

Notes

The realized measure provides information beyond what the GARCH filter can extract from squared returns alone, leading to improved volatility forecasts – especially at shorter horizons.

Reference: Hansen, P.R., Huang, Z., and Shek, H.H. (2012). “Realized GARCH: A Joint Model for Returns and Realized Measures of Volatility.” Journal of Applied Econometrics, 27(6), 877–906.

See also

garch_fit: Standard GARCH without realized measures. garch_forecast: Multi-step forecasting.

news_impact_curve(returns, model_type='GARCH', p=1, q=1, n_points=100, shock_range=3.0, dist='normal')[source]

Compute the news impact curve for a fitted GARCH-family model.

The news impact curve shows how the next-period conditional variance responds to shocks of different sizes and signs. For symmetric models (GARCH), the curve is a parabola centered at zero. For asymmetric models (EGARCH, GJR-GARCH), the curve is steeper for negative shocks, visualizing the leverage effect.

Parameters:

  • returns (Series | ndarray) – Return series for model estimation.

  • model_type (str, default: 'GARCH') –

    Type of GARCH model. Options:

    • "GARCH": Standard symmetric model.

    • "EGARCH": Nelson’s exponential GARCH.

    • "GJR": GJR-GARCH (Glosten-Jagannathan-Runkle).

  • p (int, default: 1) – GARCH lag order.

  • q (int, default: 1) – ARCH lag order.

  • n_points (int, default: 100) – Number of shock values to evaluate. Default 100.

  • shock_range (float, default: 3.0) – Range of shocks in units of unconditional standard deviation. Default 3.0 (evaluates from -3*sigma to +3*sigma).

  • dist (str, default: 'normal') – Error distribution for the model fit.

Returns:

  • shocks (np.ndarray) – Array of shock values (epsilon) evaluated, shape (n_points,).

  • conditional_variance (np.ndarray) – Corresponding next-period conditional variance for each shock.

  • model_type (str) – Type of model fitted.

  • params (dict) – Fitted model parameters.

Return type:

dict[str, Any]

Example

>>> import numpy as np
>>> from wraquant.vol.models import news_impact_curve
>>> rng = np.random.default_rng(42)
>>> returns = rng.normal(0, 1, 1000)
>>> nic = news_impact_curve(returns, model_type="GJR")
>>> nic['shocks'].shape == nic['conditional_variance'].shape
True

Notes

The news impact curve was introduced by Engle and Ng (1993) as a diagnostic tool for comparing the asymmetric response of different volatility models to news (shocks).

For GARCH(1,1): \(\sigma_{t+1}^2 = \omega + \alpha \epsilon_t^2 + \beta \sigma_t^2\) (symmetric parabola)

For GJR-GARCH: the curve has a steeper slope for \(\epsilon_t < 0\) when \(\gamma > 0\).

Reference: Engle, R.F. and Ng, V.K. (1993). “Measuring and Testing the Impact of News on Volatility.” Journal of Finance, 48(5), 1749–1778.

See also

garch_fit: Symmetric GARCH estimation. egarch_fit: Asymmetric EGARCH estimation. gjr_garch_fit: GJR-GARCH estimation.

volatility_persistence(params=None, *, alpha=None, beta=None, gamma=None, omega=None)[source]

Compute volatility persistence metrics from GARCH parameters.

Computes the half-life of volatility shocks, the persistence parameter, and the unconditional (long-run) variance implied by the GARCH model.

Parameters:

  • params (dict[str, Any] | None, default: None) – Dictionary of fitted GARCH parameters (e.g., from garch_fit()). If provided, alpha, beta, gamma, and omega are extracted automatically.

  • alpha (float | None, default: None) – Sum of ARCH coefficients. Overrides params if both are provided.

  • beta (float | None, default: None) – Sum of GARCH coefficients. Overrides params if both are provided.

  • gamma (float | None, default: None) – Sum of asymmetric (GJR) coefficients. Default 0.

  • omega (float | None, default: None) – GARCH intercept. Required for unconditional variance.

Returns:

  • persistence (float) – alpha + beta + 0.5*gamma. Values near 1 imply very persistent volatility.

  • half_life (float) – Periods for a shock to decay by 50%. ln(0.5) / ln(persistence).

  • unconditional_variance (float) – Long-run variance omega / (1 - persistence). Infinity if persistence >= 1.

  • unconditional_volatility (float) – Square root of unconditional variance.

  • mean_reversion_speed (float) – 1 - persistence. Higher values mean faster reversion.

Return type:

dict[str, float]

Example

>>> from wraquant.vol.models import volatility_persistence
>>> vp = volatility_persistence(alpha=0.05, beta=0.93, omega=0.01)
>>> vp['persistence']
0.98
>>> vp['half_life'] > 30
True

Notes

For standard GARCH, persistence = alpha + beta. For GJR-GARCH, persistence = alpha + beta + 0.5 * gamma. For EGARCH, persistence is the sum of beta coefficients.

A persistence > 1 implies non-stationarity (IGARCH behavior). Most equity series have persistence in [0.95, 0.999].

See also

garch_fit: Estimate GARCH parameters. garch_forecast: Use persistence for multi-step forecasting.

hawkes_process(events, *, max_time=None, n_iter=500)[source]

Fit a univariate Hawkes (self-exciting) point process via MLE.

The Hawkes process models volatility clustering as a self-exciting point process where each event temporarily increases the intensity (rate) of future events. This is directly analogous to how large price moves tend to cluster in financial markets.

The intensity (conditional rate) is:

\[\lambda(t) = \mu + \sum_{t_i < t} \alpha \exp(-\beta (t - t_i))\]

where \(\mu\) is the background intensity, \(\alpha\) is the jump size per event, and \(\beta\) is the exponential decay rate.

Parameters:

  • events (ndarray) – Array of event times (e.g., times of large returns or trades). Must be sorted in ascending order.

  • max_time (float | None, default: None) – End of observation window. If None, uses max(events) * 1.01.

  • n_iter (int, default: 500) – Maximum iterations for the optimizer. Default 500.

Returns:

  • mu (float) – Background (baseline) intensity.

  • alpha (float) – Excitation magnitude. Larger alpha means each event causes a bigger spike in future event intensity.

  • beta (float) – Decay rate. Larger beta means faster return to baseline intensity after an event.

  • branching_ratio (float) – alpha / beta. Must be < 1 for the process to be stationary. Values near 1 indicate strong self-excitation (volatility clustering).

  • half_life (float) – Time for excitation to decay by 50%: ln(2) / beta.

  • log_likelihood (float) – Maximized log-likelihood.

  • intensity (np.ndarray) – Estimated intensity function evaluated at each event time.

Return type:

dict[str, Any]

Example

>>> import numpy as np
>>> from wraquant.vol.models import hawkes_process
>>> rng = np.random.default_rng(42)
>>> # Simulate clustered events
>>> events = np.sort(rng.exponential(1, 100).cumsum())
>>> result = hawkes_process(events)
>>> result['mu'] > 0
True

Notes

The branching ratio alpha/beta controls the degree of self-excitation. A ratio near 1 means the process is barely stationary and events strongly cluster. In financial markets, branching ratios of 0.5-0.8 are common for trade arrivals.

The Hawkes process is equivalent to an inhomogeneous Poisson process whose intensity depends on the history of past events. It was originally developed for earthquake modeling (Hawkes, 1971) and has been widely applied to high-frequency finance.

Reference: Hawkes, A.G. (1971). “Spectra of some self-exciting and mutually exciting point processes.” Biometrika, 58(1), 83–90.

See also

stochastic_vol_sv: Continuous-time stochastic volatility model. gaussian_mixture_vol: Regime-based volatility clustering.

stochastic_vol_sv(returns, n_particles=500, n_iter=20)[source]

Fit a basic stochastic volatility model via particle filter.

The discrete-time stochastic volatility (SV) model treats log-volatility as a latent AR(1) process:

\[\begin{split}y_t &= \exp(h_t / 2) \, \epsilon_t, \quad \epsilon_t \sim N(0, 1) \\ h_t &= \mu + \phi (h_{t-1} - \mu) + \sigma_\eta \eta_t, \quad \eta_t \sim N(0, 1)\end{split}\]

where \(h_t\) is the log-variance, \(\phi\) controls persistence, \(\mu\) is the long-run mean of log-variance, and \(\sigma_\eta\) is the volatility of volatility.

Parameters:

  • returns (Series | ndarray) – Return series (not prices).

  • n_particles (int, default: 500) – Number of particles for the bootstrap particle filter. More particles = more accurate but slower. Default 500.

  • n_iter (int, default: 20) – Number of iterations for parameter learning via particle marginal MH or iterative filtering. Default 20.

Returns:

  • mu (float) – Long-run mean of log-variance.

  • phi (float) – Persistence of log-variance AR(1). Values near 1 indicate very persistent volatility.

  • sigma_eta (float) – Volatility of volatility.

  • filtered_volatility (np.ndarray) – Filtered estimate of exp(h_t/2), the conditional standard deviation.

  • log_variance (np.ndarray) – Filtered log-variance h_t.

  • log_likelihood (float) – Approximate log-likelihood from the particle filter.

Return type:

dict[str, Any]

Example

>>> import numpy as np
>>> from wraquant.vol.models import stochastic_vol_sv
>>> rng = np.random.default_rng(42)
>>> returns = rng.normal(0, 0.01, 500)
>>> result = stochastic_vol_sv(returns, n_particles=200, n_iter=5)
>>> result['phi'] > 0
True

Notes

The SV model is an alternative to GARCH that treats volatility as an unobserved latent process rather than a deterministic function of past returns. This is arguably more realistic but harder to estimate.

This implementation uses a bootstrap particle filter (sequential Monte Carlo) for likelihood evaluation and a simple grid-based parameter search. For production Bayesian SV estimation, consider PyMC or Stan.

Reference: Taylor, S.J. (1982). “Financial Returns Modelled by the Product of Two Stochastic Processes.” In Time Series Analysis: Theory and Practice 1, O.D. Anderson (ed.), 203–226.

See also

garch_fit: Deterministic volatility filtering (GARCH). gaussian_mixture_vol: Regime-switching volatility. hawkes_process: Point process for event clustering.

gaussian_mixture_vol(returns, n_components=2, random_state=42)[source]

Fit a Gaussian Mixture Model for regime-dependent volatility.

Models the return distribution as a mixture of Gaussian components, each representing a different volatility regime. This captures the empirical observation that returns alternate between calm (low-vol) and turbulent (high-vol) periods.

\[f(r_t) = \sum_{k=1}^{K} \pi_k \, \mathcal{N}(r_t; \mu_k, \sigma_k^2)\]
Parameters:

  • returns (Series | ndarray) – Return series (not prices).

  • n_components (int, default: 2) – Number of Gaussian components (regimes). Default 2 (low-vol and high-vol). Use 3 for an additional crisis regime.

  • random_state (int, default: 42) – Random seed for reproducibility.

Returns:

  • weights (np.ndarray) – Mixing weights (pi_k) for each component, shape (n_components,).

  • means (np.ndarray) – Mean return in each regime.

  • volatilities (np.ndarray) – Volatility (std dev) in each regime, sorted from lowest to highest.

  • regime_probabilities (np.ndarray) – Posterior probability of each regime for each observation, shape (T, n_components).

  • regime_labels (np.ndarray) – Most likely regime for each observation, shape (T,).

  • aic (float) – Akaike Information Criterion.

  • bic (float) – Bayesian Information Criterion.

  • model – Fitted sklearn GaussianMixture object.

Return type:

dict[str, Any]

Example

>>> import numpy as np
>>> from wraquant.vol.models import gaussian_mixture_vol
>>> rng = np.random.default_rng(42)
>>> low_vol = rng.normal(0, 0.005, 300)
>>> high_vol = rng.normal(0, 0.02, 200)
>>> returns = np.concatenate([low_vol, high_vol])
>>> result = gaussian_mixture_vol(returns, n_components=2)
>>> len(result['volatilities']) == 2
True

Notes

The GMM approach is simpler than hidden Markov models but does not model transition dynamics between regimes. For Markov regime-switching, see wraquant.regimes.

The number of components can be selected via AIC/BIC by fitting models with different n_components values.

See also

stochastic_vol_sv: Continuous latent volatility process. garch_fit: GARCH-based volatility modeling.

vol_surface_svi(strikes, forward, total_implied_var, time_to_expiry)[source]

Fit the SVI parameterization to an implied volatility smile.

The Stochastic Volatility Inspired (SVI) model of Gatheral (2004) parameterizes the total implied variance as a function of log-moneyness:

\[w(k) = a + b \left( \rho (k - m) + \sqrt{(k - m)^2 + \sigma^2} \right)\]

where \(k = \ln(K/F)\) is log-moneyness, and the five parameters (a, b, rho, m, sigma) control the level, slope, curvature, translation, and minimum variance of the smile.

Parameters:

  • strikes (ndarray) – Array of strike prices.

  • forward (float) – Forward price of the underlying.

  • total_implied_var (ndarray) – Array of total implied variances (IV^2 * T) corresponding to each strike. Same length as strikes.

  • time_to_expiry (float) – Time to expiry in years.

Returns:

  • params (dict) – Fitted SVI parameters: a (level), b (slope), rho (rotation, -1 to 1), m (translation), sigma (smoothing, > 0).

  • fitted_total_var (np.ndarray) – Fitted total implied variance at each strike.

  • fitted_iv (np.ndarray) – Fitted implied volatility at each strike (sqrt(w / T)).

  • residuals (np.ndarray) – Fitting residuals.

  • rmse (float) – Root mean squared error of the fit.

Return type:

dict[str, Any]

Example

>>> import numpy as np
>>> from wraquant.vol.models import vol_surface_svi
>>> strikes = np.array([90, 95, 100, 105, 110], dtype=float)
>>> forward = 100.0
>>> iv = np.array([0.25, 0.22, 0.20, 0.22, 0.25])
>>> T = 0.25
>>> total_var = iv**2 * T
>>> result = vol_surface_svi(strikes, forward, total_var, T)
>>> result['rmse'] < 0.01
True

Notes

The SVI parameterization is widely used by practitioners for interpolating and extrapolating the volatility smile. It naturally produces realistic smile shapes and can be calibrated to satisfy no-arbitrage constraints (Gatheral and Jacquier, 2014).

Key properties:

  • rho < 0: negative skew (typical for equities).

  • b controls the overall slope of the wings.

  • sigma controls the curvature near the money.

Reference: Gatheral, J. (2004). “A Parsimonious Arbitrage-Free Implied Volatility Parameterization.” Presentation at Global Derivatives & Risk Management, Madrid.

See also

variance_risk_premium: Compute VRP from implied and realized vol.

variance_risk_premium(implied_vol, realized_vol, annualized=True)[source]

Compute the variance risk premium (VRP).

The variance risk premium is the difference between risk-neutral (implied) and physical (realized) variance. A positive VRP indicates that investors pay a premium for volatility protection, which is the empirical norm for equity indices.

\[VRP_t = IV_t^2 - RV_t^2\]

A consistently positive VRP is the economic basis for volatility selling strategies (e.g., selling straddles, variance swaps).

Parameters:

  • implied_vol (Series | ndarray) – Implied volatility series (e.g., VIX / 100 for S&P 500). Should be in the same units as realized_vol (both annualized or both not).

  • realized_vol (Series | ndarray) – Realized volatility series. Same length and alignment as implied_vol.

  • annualized (bool, default: True) – Whether the input volatilities are annualized. If True, the VRP is in annualized variance units.

Returns:

  • vrp (np.ndarray) – Variance risk premium time series. Positive values = implied > realized (normal).

  • mean_vrp (float) – Average VRP over the sample.

  • vrp_ratio (np.ndarray) – IV / RV ratio. Values > 1 indicate implied exceeds realized.

  • pct_positive (float) – Fraction of observations where VRP > 0. Typically 60–80% for equity indices.

  • vol_spread (np.ndarray) – Simple difference IV - RV (in volatility, not variance units).

Return type:

dict[str, Any]

Example

>>> import numpy as np
>>> from wraquant.vol.models import variance_risk_premium
>>> iv = np.array([0.20, 0.22, 0.18, 0.25, 0.19])
>>> rv = np.array([0.15, 0.17, 0.16, 0.20, 0.14])
>>> result = variance_risk_premium(iv, rv)
>>> result['mean_vrp'] > 0
True

Notes

The VRP is one of the most robust risk premiums in financial markets. It reflects the cost of insurance against volatility spikes and is related to jump risk, model uncertainty, and aggregate risk aversion.

Reference: Bollerslev, T., Tauchen, G., and Zhou, H. (2009). “Expected Stock Returns and Variance Risk Premia.” Review of Financial Studies, 22(11), 4463–4492.

See also

garch_fit: Estimate realized conditional volatility. vol_surface_svi: Implied volatility surface fitting.

GARCH Models

Volatility models — GARCH family, stochastic vol, and Hawkes processes.

This module provides deep integration with the arch library for GARCH-family models, plus pure implementations of stochastic volatility and Hawkes processes for modeling volatility clustering and self-exciting dynamics.

All GARCH-family functions return rich result dictionaries including fitted parameters, conditional volatility series, diagnostic statistics, and the underlying model object for further analysis.

ewma_volatility(returns, span=30, annualize=True, periods_per_year=252)[source]

Exponentially weighted moving average volatility.

The EWMA model of J.P. Morgan’s RiskMetrics (1996) assigns exponentially decaying weights to past squared returns. The decay factor lambda = 1 - 2/(span+1).

\[\sigma_t^2 = \lambda \sigma_{t-1}^2 + (1 - \lambda) r_{t-1}^2\]
Parameters:

  • returns (Series) – Return series (not prices).

  • span (int, default: 30) – EWMA span parameter. The decay factor is computed as lambda = 1 - 2/(span+1). RiskMetrics uses span ~= 73 (lambda = 0.94 for daily data). Default 30.

  • annualize (bool, default: True) – Whether to annualize the volatility by multiplying by sqrt(periods_per_year).

  • periods_per_year (int, default: 252) – Trading periods per year. Use 252 for daily, 52 for weekly, 12 for monthly.

Return type:

Series

Returns:

EWMA volatility series with same index as input.

Example

>>> import pandas as pd, numpy as np
>>> rng = np.random.default_rng(42)
>>> returns = pd.Series(rng.normal(0, 0.01, 500))
>>> vol = ewma_volatility(returns, span=30)
>>> vol.iloc[-1] > 0
True

Notes

RiskMetrics (1996) recommended lambda = 0.94 for daily data and lambda = 0.97 for monthly data. The EWMA model is a restricted IGARCH(1,1) with zero intercept.

See also

garch_fit: More flexible parametric volatility model. realized_volatility: Non-parametric rolling window estimator.

garch_fit(returns, p=1, q=1, mean='Constant', dist='normal', **kwargs)[source]

Fit a GARCH(p, q) model to a return series.

The Generalized Autoregressive Conditional Heteroskedasticity model of Bollerslev (1986) captures volatility clustering – the tendency for large moves to follow large moves. The conditional variance is:

\[\sigma_t^2 = \omega + \sum_{i=1}^{p} \alpha_i \epsilon_{t-i}^2 + \sum_{j=1}^{q} \beta_j \sigma_{t-j}^2\]
Parameters:

  • returns (Series | ndarray) – Return series (not prices). Should be percentage returns or log returns, typically in the range [-10, 10].

  • p (int, default: 1) – Order of the GARCH (lagged variance) terms. Higher values capture longer volatility memory. Default 1 is standard.

  • q (int, default: 1) – Order of the ARCH (lagged squared residual) terms. Default 1.

  • mean (str, default: 'Constant') –

    Mean model specification. Options:

    • "Constant": constant mean (most common)

    • "Zero": zero mean (for demeaned returns)

    • "ARX": autoregressive with exogenous variables

  • dist (str, default: 'normal') –

    Error distribution. Options:

    • "normal": Gaussian (fastest, least flexible)

    • "t": Student’s t (captures fat tails)

    • "skewt": Hansen’s skewed t (tails + asymmetry)

    • "ged": Generalized Error Distribution

  • **kwargs (Any) – Additional keyword arguments passed to arch.arch_model() (e.g., lags for ARX mean).

Returns:

  • model_name (str) – Model identifier, e.g. "GARCH(1,1)".

  • params (dict) – Fitted parameters (omega, alpha, beta, plus distribution params if non-normal).

  • conditional_volatility (pd.Series) – Time series of estimated conditional standard deviation. Multiply by sqrt(252) to annualize if returns are daily.

  • standardized_residuals (np.ndarray) – Residuals divided by conditional std dev. Should be ~i.i.d. if model is correct.

  • aic (float) – Akaike Information Criterion (lower = better).

  • bic (float) – Bayesian Information Criterion.

  • log_likelihood (float) – Maximized log-likelihood.

  • persistence (float) – Sum of alpha + beta. Values near 1 indicate high volatility persistence (IGARCH if exactly 1).

  • half_life (float) – Number of periods for a volatility shock to decay to half its initial impact.

  • unconditional_variance (float) – Long-run variance omega / (1 - alpha - beta).

  • ljung_box (dict) – Ljung-Box test on squared standardized residuals. If p-value > 0.05, model captures vol clustering.

  • model – The fitted arch model result object for further analysis (forecasting, simulation, etc.).

Return type:

dict[str, Any]

Example

>>> import numpy as np
>>> from wraquant.vol.models import garch_fit
>>> rng = np.random.default_rng(42)
>>> returns = rng.normal(0, 1, 1000)
>>> result = garch_fit(returns, p=1, q=1, dist="normal")
>>> 0 < result['persistence'] < 2
True

Notes

The GARCH(1,1) model is the workhorse of volatility modeling. Persistence (alpha + beta) near 1 is common for financial returns. If persistence > 1, consider IGARCH or FIGARCH.

Returns are internally scaled by 100 before fitting to improve numerical stability (the arch library convention). All outputs are rescaled back to the original units.

Reference: Bollerslev, T. (1986). “Generalized Autoregressive Conditional Heteroskedasticity.” Journal of Econometrics, 31, 307–327.

See also

egarch_fit: Asymmetric GARCH capturing leverage effect. gjr_garch_fit: Alternative asymmetric specification. garch_forecast: Multi-step ahead volatility forecasting. news_impact_curve: Visualize asymmetric shock response.

egarch_fit(returns, p=1, q=1, mean='Constant', dist='normal', **kwargs)[source]

Fit an EGARCH(p, q) model (exponential GARCH).

The EGARCH model of Nelson (1991) captures the leverage effect – negative shocks increase volatility more than positive shocks of equal magnitude. The model specifies the log of conditional variance, so positivity of variance is guaranteed without parameter constraints:

\[\ln \sigma_t^2 = \omega + \sum_{i=1}^{q} \alpha_i \left( \left|z_{t-i}\right| - E\left|z_{t-i}\right| \right) + \sum_{i=1}^{q} \gamma_i z_{t-i} + \sum_{j=1}^{p} \beta_j \ln \sigma_{t-j}^2\]

where \(z_t = \epsilon_t / \sigma_t\) are standardized residuals and \(\gamma_i\) captures the asymmetric (leverage) response.

Parameters:

  • returns (Series | ndarray) – Return series (not prices).

  • p (int, default: 1) – Order of GARCH (lagged log-variance) terms.

  • q (int, default: 1) – Order of ARCH (asymmetric innovation) terms.

  • mean (str, default: 'Constant') – Mean model specification ("Constant", "Zero", "ARX").

  • dist (str, default: 'normal') – Error distribution ("normal", "t", "skewt", "ged").

  • **kwargs (Any) – Additional keyword arguments passed to arch.arch_model().

Return type:

dict[str, Any]

Returns:

Dictionary with the same structure as garch_fit(). The params dict includes gamma terms capturing the leverage effect. A negative gamma means negative shocks increase volatility more than positive shocks.

Example

>>> import numpy as np
>>> from wraquant.vol.models import egarch_fit
>>> rng = np.random.default_rng(42)
>>> returns = rng.normal(0, 1, 1000)
>>> result = egarch_fit(returns, p=1, q=1)
>>> result['conditional_volatility'].iloc[-1] > 0
True

Notes

The EGARCH model is preferred when:

  1. Leverage effects are important (equity markets).

  2. Parameter positivity constraints are problematic.

The log specification means persistence is measured by the sum of beta coefficients. Unconditional variance may not have a closed form depending on the distribution.

Reference: Nelson, D.B. (1991). “Conditional Heteroskedasticity in Asset Returns: A New Approach.” Econometrica, 59(2), 347–370.

See also

garch_fit: Symmetric GARCH model. gjr_garch_fit: Alternative asymmetric GARCH. news_impact_curve: Compare asymmetric responses across models.

gjr_garch_fit(returns, p=1, q=1, o=1, mean='Constant', dist='normal', **kwargs)[source]

Fit a GJR-GARCH(p, o, q) model (Glosten-Jagannathan-Runkle).

The GJR-GARCH model of Glosten, Jagannathan, and Runkle (1993) augments the standard GARCH with an indicator function that activates on negative shocks, capturing the leverage effect:

\[\sigma_t^2 = \omega + \sum_{i=1}^{q} \alpha_i \epsilon_{t-i}^2 + \sum_{i=1}^{o} \gamma_i \epsilon_{t-i}^2 I(\epsilon_{t-i} < 0) + \sum_{j=1}^{p} \beta_j \sigma_{t-j}^2\]

where \(I(\cdot)\) is the indicator function. The gamma coefficient captures the extra volatility increase from negative shocks (bad news).

Parameters:

  • returns (Series | ndarray) – Return series (not prices).

  • p (int, default: 1) – Order of GARCH (lagged variance) terms.

  • q (int, default: 1) – Order of ARCH (lagged squared residual) terms.

  • o (int, default: 1) – Order of asymmetric (threshold) terms. Default 1.

  • mean (str, default: 'Constant') – Mean model specification.

  • dist (str, default: 'normal') – Error distribution.

  • **kwargs (Any) – Additional keyword arguments passed to arch.arch_model().

Return type:

dict[str, Any]

Returns:

Dictionary with the same structure as garch_fit(). The params dict includes gamma[1] capturing the asymmetric response. Positive gamma means negative shocks increase volatility more.

The persistence for GJR-GARCH is alpha + beta + 0.5 * gamma (assuming symmetric distribution).

Example

>>> import numpy as np
>>> from wraquant.vol.models import gjr_garch_fit
>>> rng = np.random.default_rng(42)
>>> returns = rng.normal(0, 1, 1000)
>>> result = gjr_garch_fit(returns, p=1, q=1)
>>> 'gamma[1]' in result['params'] or 'gamma' in str(result['params'])
True

Notes

The GJR-GARCH is equivalent to TARCH (Threshold ARCH) when applied to the variance (not standard deviation). It is the default asymmetric model in many risk management applications.

For stationarity, the persistence alpha + beta + 0.5 * gamma < 1 is required.

Reference: Glosten, L.R., Jagannathan, R., and Runkle, D.E. (1993). “On the Relation between the Expected Value and the Volatility of the Nominal Excess Return on Stocks.” Journal of Finance, 48(5), 1779–1801.

See also

egarch_fit: Log-variance based asymmetric model. garch_fit: Symmetric baseline GARCH. news_impact_curve: Compare asymmetric responses.

figarch_fit(returns, p=1, q=1, mean='Constant', dist='normal', **kwargs)[source]

Fit a FIGARCH(p, d, q) model (Fractionally Integrated GARCH).

The FIGARCH model of Baillie, Bollerslev, and Mikkelsen (1996) allows for long memory in the conditional variance process. The fractional differencing parameter d is estimated from the data and captures the slow hyperbolic decay of volatility autocorrelations that is frequently observed in financial time series.

\[\sigma_t^2 = \omega + \beta \sigma_{t-1}^2 + [1 - \beta L - (1-L)^d (1 - \alpha L)] \epsilon_t^2\]

where d in (0, 1) is the fractional integration parameter, and L is the lag operator.

Parameters:

  • returns (Series | ndarray) – Return series (not prices).

  • p (int, default: 1) – Order of GARCH terms (typically 1).

  • q (int, default: 1) – Order of ARCH terms (typically 1).

  • mean (str, default: 'Constant') – Mean model specification.

  • dist (str, default: 'normal') – Error distribution.

  • **kwargs (Any) – Additional keyword arguments passed to arch.arch_model().

Return type:

dict[str, Any]

Returns:

Dictionary with the same structure as garch_fit(). The params dict includes the fractional integration parameter d. Values of d near 0 imply short memory (like GARCH), values near 0.5 imply long memory.

Example

>>> import numpy as np
>>> from wraquant.vol.models import figarch_fit
>>> rng = np.random.default_rng(42)
>>> returns = rng.normal(0, 1, 1000)
>>> result = figarch_fit(returns, p=1, q=1)
>>> result['aic'] < 0 or result['aic'] > 0  # AIC is finite
True

Notes

FIGARCH nests GARCH (d=0) and IGARCH (d=1). Empirically, most financial return series have d in (0.3, 0.5), suggesting long memory in volatility.

The arch library implements FIGARCH via the "FIGARCH" volatility specification. The power parameter defaults to 2.0 (standard FIGARCH).

Reference: Baillie, R.T., Bollerslev, T., and Mikkelsen, H.O. (1996). “Fractionally Integrated Generalized Autoregressive Conditional Heteroskedasticity.” Journal of Econometrics, 74(1), 3–30.

See also

garch_fit: Standard GARCH without long memory. harch_fit: Heterogeneous ARCH as alternative long-memory model. volatility_persistence: Measure persistence and half-life.

aparch_fit(returns, p=1, o=1, q=1, mean='Constant', dist='normal', **kwargs)[source]

Fit an APARCH(p, o, q) model (Asymmetric Power ARCH).

The APARCH model of Ding, Granger, and Engle (1993) generalises the standard GARCH by raising the conditional standard deviation to a freely estimated power delta and introducing a leverage (rotation) parameter gamma for each asymmetric term:

\[\sigma_t^\delta = \omega + \sum_{i=1}^{q} \alpha_i \left( |\epsilon_{t-i}| - \gamma_i \epsilon_{t-i} \right)^\delta + \sum_{j=1}^{p} \beta_j \sigma_{t-j}^\delta\]

The power parameter delta is estimated from the data rather than being fixed at 2 (GARCH) or 1 (TARCH). This gives the model extra flexibility to capture the shape of the conditional distribution of volatility. When delta = 2 the model collapses to GJR-GARCH; when delta = 1 it collapses to TARCH.

Parameters:

  • returns (Series | ndarray) – Return series (not prices). Should be percentage returns or log returns.

  • p (int, default: 1) – Order of GARCH (lagged powered standard deviation) terms.

  • o (int, default: 1) – Order of asymmetric (leverage) terms. Default 1.

  • q (int, default: 1) – Order of ARCH (lagged powered absolute residual) terms.

  • mean (str, default: 'Constant') –

    Mean model specification. Options:

    • "Constant": constant mean (most common)

    • "Zero": zero mean (for demeaned returns)

    • "ARX": autoregressive with exogenous variables

  • dist (str, default: 'normal') –

    Error distribution. Options:

    • "normal": Gaussian

    • "t": Student’s t (fat tails)

    • "skewt": Hansen’s skewed t

    • "ged": Generalized Error Distribution

  • **kwargs (Any) – Additional keyword arguments passed to arch.arch_model().

Returns:

  • model_name (str) – "APARCH(p,o,q)".

  • params (dict) – Fitted parameters including omega, alpha[i], gamma[i] (leverage), beta[j], and the estimated power parameter delta (often labelled delta or embedded in the volatility process params).

  • conditional_volatility (pd.Series) – Conditional standard deviation series.

  • standardized_residuals (np.ndarray) – Model residuals divided by conditional std dev.

  • aic / bic / log_likelihood (float) – Fit quality metrics.

  • persistence (float) – Volatility persistence.

  • half_life (float) – Half-life of volatility shocks.

  • unconditional_variance (float) – Long-run variance.

  • ljung_box (dict) – Ljung-Box test on squared standardized residuals.

  • model – Fitted arch result object.

Return type:

dict[str, Any]

Example

>>> import numpy as np
>>> from wraquant.vol.models import aparch_fit
>>> rng = np.random.default_rng(42)
>>> returns = rng.normal(0, 1, 1000)
>>> result = aparch_fit(returns, p=1, o=1, q=1)
>>> result['conditional_volatility'].iloc[-1] > 0
True

Notes

The key advantage of APARCH is the estimated power parameter delta. Empirical studies (Ding et al. 1993, Laurent 2004) find that delta is often close to 1 for daily equity returns rather than 2, implying the conditional standard deviation (not variance) follows a more parsimonious dynamic.

APARCH nests at least seven other ARCH-class models as special cases, making it useful for model selection via the estimated delta.

Reference: Ding, Z., Granger, C.W.J., and Engle, R.F. (1993). “A Long Memory Property of Stock Market Returns and a New Model.” Journal of Empirical Finance, 1, 83–106.

See also

garch_fit: Standard GARCH (delta fixed at 2, no leverage). gjr_garch_fit: GJR-GARCH (delta fixed at 2, with leverage). egarch_fit: Log-variance based asymmetric model.

garch_forecast(returns, p=1, q=1, mean='Constant', dist='normal', horizon=10, method='analytic', simulations=1000, **kwargs)[source]

Multi-step ahead volatility forecast from a GARCH model.

Fits a GARCH(p,q) model and produces multi-step forecasts of conditional variance, with optional confidence intervals via simulation.

Parameters:

  • returns (Series | ndarray) – Return series (not prices).

  • p (int, default: 1) – GARCH lag order.

  • q (int, default: 1) – ARCH lag order.

  • mean (str, default: 'Constant') – Mean model specification.

  • dist (str, default: 'normal') – Error distribution.

  • horizon (int, default: 10) – Number of steps ahead to forecast. Default 10.

  • method (str, default: 'analytic') –

    Forecasting method. Options:

    • "analytic": Closed-form multi-step variance forecast. Fastest and exact for GARCH. Does not provide intervals.

    • "simulation": Monte Carlo simulation of future paths. Provides confidence intervals but slower.

  • simulations (int, default: 1000) – Number of simulation paths (only used when method="simulation"). Default 1000.

  • **kwargs (Any) – Additional keyword arguments passed to arch.arch_model().

Returns:

  • model_name (str) – "GARCH(p,q) forecast".

  • forecast_variance (np.ndarray) – Forecasted conditional variance for each horizon step, shape (horizon,).

  • forecast_volatility (np.ndarray) – Square root of forecast variance, shape (horizon,).

  • confidence_intervals (dict | None) – If method is "simulation", contains "lower_5" and "upper_95" arrays. None for analytic method.

  • fit_result (dict) – Full GARCH fit result from garch_fit().

Return type:

dict[str, Any]

Example

>>> import numpy as np
>>> from wraquant.vol.models import garch_forecast
>>> rng = np.random.default_rng(42)
>>> returns = rng.normal(0, 1, 1000)
>>> fc = garch_forecast(returns, horizon=5)
>>> len(fc['forecast_volatility']) == 5
True

Notes

Analytic forecasts assume the model is correctly specified and use the recursive formula for multi-step variance:

\[E[\sigma_{T+h}^2 | \mathcal{F}_T] = \omega \frac{1 - (\alpha + \beta)^h}{1 - \alpha - \beta} + (\alpha + \beta)^h \sigma_{T+1}^2\]

Simulation forecasts draw from the fitted error distribution and roll the GARCH recursion forward, providing a distribution of future paths.

See also

garch_fit: Fit the underlying GARCH model. volatility_persistence: Analyze shock decay properties.

garch_rolling_forecast(returns, window=None, horizon=1, vol='GARCH', p=1, q=1, dist='normal', refit_every=1, **kwargs)[source]

Rolling or expanding window 1-step-ahead GARCH volatility forecast.

At each step the model is fit on past data (fixed window or expanding) and a horizon-step-ahead forecast is produced. The resulting series of out-of-sample forecasts can be compared against realised volatility for back-testing purposes.

\[\hat{\sigma}_{t+1|t}^2 = \omega + \alpha \epsilon_t^2 + \beta \sigma_t^2\]

The forecast is produced at each date t using only data up to t, providing a genuinely out-of-sample volatility estimate.

Parameters:

  • returns (Series | ndarray) – Return series. If a pd.Series with a DatetimeIndex, dates are preserved in the output.

  • window (int | None, default: None) – Rolling window size (number of observations). If None an expanding window is used, starting from the first 100 observations. A fixed window is common in practice (e.g., 500 or 1000 days).

  • horizon (int, default: 1) – Forecast horizon (steps ahead). Default 1 for 1-step.

  • vol (str, default: 'GARCH') – Volatility model type. Any valid arch vol spec: "GARCH", "EGARCH", "GJR-GARCH", "APARCH", etc. Default "GARCH".

  • p (int, default: 1) – GARCH order.

  • q (int, default: 1) – ARCH order.

  • dist (str, default: 'normal') – Error distribution ("normal", "t", "skewt").

  • refit_every (int, default: 1) – Refit the model every n steps. Setting this > 1 greatly speeds up the rolling forecast by re-using the previous fit for intermediate steps. Default 1 (refit every step).

  • **kwargs (Any) – Additional keyword arguments passed to arch.arch_model().

Returns:

  • forecasts (pd.Series) – Series of forecasted conditional volatility (standard deviation), indexed by the date at which the forecast was made.

  • actuals (pd.Series) – Absolute returns as a proxy for realised volatility, aligned with forecasts.

  • dates (pd.Index | np.ndarray) – Forecast dates or integer positions.

Return type:

dict[str, Any]

Example

>>> import numpy as np, pandas as pd
>>> from wraquant.vol.models import garch_rolling_forecast
>>> rng = np.random.default_rng(42)
>>> ret = pd.Series(rng.normal(0, 0.01, 300))
>>> result = garch_rolling_forecast(ret, window=200, refit_every=10)
>>> len(result['forecasts']) > 0
True
>>> (result['forecasts'] > 0).all()
True

Notes

Setting refit_every > 1 avoids refitting the model at every single step. Between refits the most recent parameter estimates are used and the conditional variance is simply rolled forward using the GARCH recursion. This is standard practice in the literature (Hansen & Lunde 2005) and can reduce computation by an order of magnitude.

When vol="EGARCH" the arch library fits the model on log variance, so internal rescaling differs. The output is always in the same units as the input returns.

Reference: Hansen, P.R. and Lunde, A. (2005). “A Forecast Comparison of Volatility Models: Does Anything Beat a GARCH(1,1)?” Journal of Applied Econometrics, 20, 873–889.

See also

garch_forecast: Single multi-step forecast from a fitted model. garch_model_selection: Compare GARCH variants systematically.

garch_model_selection(returns, p=1, q=1, distributions=None, mean='Constant', **kwargs)[source]

Compare GARCH, EGARCH, and GJR-GARCH across error distributions.

Fits every combination of volatility specification and error distribution, then ranks the models by AIC (and reports BIC, log-likelihood, and persistence for each).

This is the standard model-selection exercise recommended by Hansen & Lunde (2005) to determine which GARCH variant and error distribution best describes a given return series.

Parameters:

  • returns (Series | ndarray) – Return series (not prices).

  • p (int, default: 1) – GARCH order (applied to all models).

  • q (int, default: 1) – ARCH order (applied to all models).

  • distributions (list[str] | None, default: None) – List of error distributions to try. Defaults to ["normal", "t", "skewt"].

  • mean (str, default: 'Constant') – Mean model specification applied to all fits.

  • **kwargs (Any) – Additional keyword arguments passed to arch.arch_model().

Returns:

  • model (str) – Model name (e.g. "EGARCH-t").

  • vol_model (str) – Volatility specification.

  • distribution (str) – Error distribution.

  • aic (float) – Akaike Information Criterion.

  • bic (float) – Bayesian Information Criterion.

  • log_likelihood (float) – Maximised log-likelihood.

  • persistence (float) – Volatility persistence.

  • num_params (int) – Number of estimated parameters.

Return type:

DataFrame

Example

>>> import numpy as np
>>> from wraquant.vol.models import garch_model_selection
>>> rng = np.random.default_rng(42)
>>> returns = rng.normal(0, 1, 1000)
>>> df = garch_model_selection(returns)
>>> 'aic' in df.columns and 'bic' in df.columns
True
>>> len(df) >= 6  # 3 models x 2+ distributions
True

Notes

Models that fail to converge are silently excluded from the results table. AIC penalises complexity less than BIC, so the best AIC and best BIC model may differ – both are reported.

Reference: Hansen, P.R. and Lunde, A. (2005). “A Forecast Comparison of Volatility Models: Does Anything Beat a GARCH(1,1)?” Journal of Applied Econometrics, 20, 873–889.

See also

garch_fit: Fit a specific GARCH model. egarch_fit: Fit EGARCH. gjr_garch_fit: Fit GJR-GARCH.

dcc_fit(returns, p=1, q=1, dist='normal')[source]

Fit a DCC-GARCH model for multivariate dynamic correlations.

The Dynamic Conditional Correlation model of Engle (2002) estimates time-varying correlations between multiple asset returns using a two-step procedure:

  1. Fit univariate GARCH(p,q) to each series.

  2. Estimate DCC parameters (a, b) via MLE on standardized residuals.

The DCC correlation dynamics are:

\[ \begin{align}\begin{aligned}Q_t = (1 - a - b) \bar{Q} + a \, z_{t-1} z_{t-1}^\top + b \, Q_{t-1}\\R_t = \text{diag}(Q_t)^{-1/2} \, Q_t \, \text{diag}(Q_t)^{-1/2}\end{aligned}\end{align} \]

where \(\bar{Q}\) is the unconditional correlation matrix of standardized residuals.

Parameters:

  • returns (DataFrame | ndarray) – Return data, either a DataFrame with k columns (one per asset) or an array of shape (T, k).

  • p (int, default: 1) – GARCH lag order for each univariate model.

  • q (int, default: 1) – ARCH lag order for each univariate model.

  • dist (str, default: 'normal') – Error distribution for univariate GARCH fits.

Returns:

  • dcc_params (dict) – DCC parameters {"a": ..., "b": ...}.

  • univariate_results (list[dict]) – Per-asset GARCH fit results from garch_fit().

  • conditional_correlations (np.ndarray) – Array of shape (T, k, k) with time-varying correlation matrices.

  • conditional_covariances (np.ndarray) – Array of shape (T, k, k) with time-varying covariance matrices.

  • conditional_volatilities (np.ndarray) – Array of shape (T, k) with per-asset conditional volatilities.

  • standardized_residuals (np.ndarray) – Array of shape (T, k).

  • qbar (np.ndarray) – Unconditional correlation matrix.

Return type:

dict[str, Any]

Example

>>> import numpy as np, pandas as pd
>>> from wraquant.vol.models import dcc_fit
>>> rng = np.random.default_rng(42)
>>> ret = pd.DataFrame({
...     'A': rng.normal(0, 0.01, 500),
...     'B': rng.normal(0, 0.015, 500),
... })
>>> result = dcc_fit(ret)
>>> result['conditional_correlations'].shape == (500, 2, 2)
True

Notes

This implementation uses the arch library for the univariate GARCH step, then applies the pure scipy/numpy DCC estimation from wraquant.risk.dcc.

The DCC model is widely used in portfolio management, value-at-risk computation, and hedging applications where time-varying correlations are important.

Reference: Engle, R. (2002). “Dynamic Conditional Correlation: A Simple Class of Multivariate Generalized Autoregressive Conditional Heteroskedasticity Models.” Journal of Business & Economic Statistics, 20(3), 339–350.

See also

garch_fit: Univariate GARCH estimation. garch_forecast: Volatility forecasting.

realized_garch(returns, realized_vol, p=1, q=1, mean='Constant', dist='normal', **kwargs)[source]

Fit a Realized GARCH model augmented with realized volatility.

The Realized GARCH model of Hansen, Huang, and Shek (2012) augments the standard GARCH with a measurement equation linking the conditional variance to a realized volatility measure, improving forecasting performance compared to standard GARCH.

This implementation uses the arch library’s ARX mean model with lagged realized volatility as an exogenous regressor, capturing the predictive information in the realized measure.

Parameters:

  • returns (Series | ndarray) – Return series (not prices).

  • realized_vol (Series | ndarray) – Corresponding realized volatility measure (must have the same length as returns). Can be any realized measure: realized variance, Parkinson, Garman-Klass, etc.

  • p (int, default: 1) – GARCH lag order.

  • q (int, default: 1) – ARCH lag order.

  • mean (str, default: 'Constant') – Mean model specification. Default "Constant" uses a constant plus the realized vol regressor.

  • dist (str, default: 'normal') – Error distribution.

  • **kwargs (Any) – Additional keyword arguments passed to arch.arch_model().

Return type:

dict[str, Any]

Returns:

Dictionary with the same structure as garch_fit(), plus:

  • realized_vol_used (np.ndarray) – The realized volatility series used in the model.

Example

>>> import numpy as np
>>> from wraquant.vol.models import realized_garch
>>> rng = np.random.default_rng(42)
>>> returns = rng.normal(0, 0.01, 500)
>>> rv = np.abs(returns).rolling(20).mean() if hasattr(returns, 'rolling') else np.convolve(np.abs(returns), np.ones(20)/20, mode='same')
>>> # result = realized_garch(returns, rv)

Notes

The realized measure provides information beyond what the GARCH filter can extract from squared returns alone, leading to improved volatility forecasts – especially at shorter horizons.

Reference: Hansen, P.R., Huang, Z., and Shek, H.H. (2012). “Realized GARCH: A Joint Model for Returns and Realized Measures of Volatility.” Journal of Applied Econometrics, 27(6), 877–906.

See also

garch_fit: Standard GARCH without realized measures. garch_forecast: Multi-step forecasting.

harch_fit(returns, lags=None, mean='Constant', dist='normal', **kwargs)[source]

Fit a HARCH model (Heterogeneous ARCH).

The HARCH model of Mueller et al. (1997) aggregates squared returns over multiple time horizons to capture the heterogeneous nature of market participants operating at different frequencies (intraday traders, daily, weekly, monthly):

\[\sigma_t^2 = \omega + \sum_{j=1}^{J} \alpha_j \left( \sum_{i=1}^{l_j} \epsilon_{t-i} \right)^2 / l_j\]

where \(l_1 < l_2 < \ldots < l_J\) are the lag lengths corresponding to different time horizons.

Parameters:

  • returns (Series | ndarray) – Return series (not prices).

  • lags (list[int] | None, default: None) – List of lag lengths representing different time horizons. Default is [1, 5, 22] corresponding to daily, weekly, and monthly horizons for daily data.

  • mean (str, default: 'Constant') – Mean model specification.

  • dist (str, default: 'normal') – Error distribution.

  • **kwargs (Any) – Additional keyword arguments passed to arch.arch_model().

Return type:

dict[str, Any]

Returns:

Dictionary with the same structure as garch_fit(). The params dict includes coefficients for each aggregation horizon.

Example

>>> import numpy as np
>>> from wraquant.vol.models import harch_fit
>>> rng = np.random.default_rng(42)
>>> returns = rng.normal(0, 1, 1000)
>>> result = harch_fit(returns, lags=[1, 5, 22])
>>> result['conditional_volatility'].iloc[-1] > 0
True

Notes

HARCH is particularly useful for modeling high-frequency data where different market participants operate on vastly different time scales. It provides a parsimonious alternative to high-order ARCH or GARCH models.

Reference: Mueller, U.A., Dacorogna, M.M., Dave, R.D., Olsen, R.B., Pictet, O.V., and von Weizsacker, J.E. (1997). “Volatilities of Different Time Resolutions.” Journal of Empirical Finance, 4(2-3), 213–239.

See also

garch_fit: Standard GARCH model. figarch_fit: Long-memory GARCH via fractional integration.

news_impact_curve(returns, model_type='GARCH', p=1, q=1, n_points=100, shock_range=3.0, dist='normal')[source]

Compute the news impact curve for a fitted GARCH-family model.

The news impact curve shows how the next-period conditional variance responds to shocks of different sizes and signs. For symmetric models (GARCH), the curve is a parabola centered at zero. For asymmetric models (EGARCH, GJR-GARCH), the curve is steeper for negative shocks, visualizing the leverage effect.

Parameters:

  • returns (Series | ndarray) – Return series for model estimation.

  • model_type (str, default: 'GARCH') –

    Type of GARCH model. Options:

    • "GARCH": Standard symmetric model.

    • "EGARCH": Nelson’s exponential GARCH.

    • "GJR": GJR-GARCH (Glosten-Jagannathan-Runkle).

  • p (int, default: 1) – GARCH lag order.

  • q (int, default: 1) – ARCH lag order.

  • n_points (int, default: 100) – Number of shock values to evaluate. Default 100.

  • shock_range (float, default: 3.0) – Range of shocks in units of unconditional standard deviation. Default 3.0 (evaluates from -3*sigma to +3*sigma).

  • dist (str, default: 'normal') – Error distribution for the model fit.

Returns:

  • shocks (np.ndarray) – Array of shock values (epsilon) evaluated, shape (n_points,).

  • conditional_variance (np.ndarray) – Corresponding next-period conditional variance for each shock.

  • model_type (str) – Type of model fitted.

  • params (dict) – Fitted model parameters.

Return type:

dict[str, Any]

Example

>>> import numpy as np
>>> from wraquant.vol.models import news_impact_curve
>>> rng = np.random.default_rng(42)
>>> returns = rng.normal(0, 1, 1000)
>>> nic = news_impact_curve(returns, model_type="GJR")
>>> nic['shocks'].shape == nic['conditional_variance'].shape
True

Notes

The news impact curve was introduced by Engle and Ng (1993) as a diagnostic tool for comparing the asymmetric response of different volatility models to news (shocks).

For GARCH(1,1): \(\sigma_{t+1}^2 = \omega + \alpha \epsilon_t^2 + \beta \sigma_t^2\) (symmetric parabola)

For GJR-GARCH: the curve has a steeper slope for \(\epsilon_t < 0\) when \(\gamma > 0\).

Reference: Engle, R.F. and Ng, V.K. (1993). “Measuring and Testing the Impact of News on Volatility.” Journal of Finance, 48(5), 1749–1778.

See also

garch_fit: Symmetric GARCH estimation. egarch_fit: Asymmetric EGARCH estimation. gjr_garch_fit: GJR-GARCH estimation.

volatility_persistence(params=None, *, alpha=None, beta=None, gamma=None, omega=None)[source]

Compute volatility persistence metrics from GARCH parameters.

Computes the half-life of volatility shocks, the persistence parameter, and the unconditional (long-run) variance implied by the GARCH model.

Parameters:

  • params (dict[str, Any] | None, default: None) – Dictionary of fitted GARCH parameters (e.g., from garch_fit()). If provided, alpha, beta, gamma, and omega are extracted automatically.

  • alpha (float | None, default: None) – Sum of ARCH coefficients. Overrides params if both are provided.

  • beta (float | None, default: None) – Sum of GARCH coefficients. Overrides params if both are provided.

  • gamma (float | None, default: None) – Sum of asymmetric (GJR) coefficients. Default 0.

  • omega (float | None, default: None) – GARCH intercept. Required for unconditional variance.

Returns:

  • persistence (float) – alpha + beta + 0.5*gamma. Values near 1 imply very persistent volatility.

  • half_life (float) – Periods for a shock to decay by 50%. ln(0.5) / ln(persistence).

  • unconditional_variance (float) – Long-run variance omega / (1 - persistence). Infinity if persistence >= 1.

  • unconditional_volatility (float) – Square root of unconditional variance.

  • mean_reversion_speed (float) – 1 - persistence. Higher values mean faster reversion.

Return type:

dict[str, float]

Example

>>> from wraquant.vol.models import volatility_persistence
>>> vp = volatility_persistence(alpha=0.05, beta=0.93, omega=0.01)
>>> vp['persistence']
0.98
>>> vp['half_life'] > 30
True

Notes

For standard GARCH, persistence = alpha + beta. For GJR-GARCH, persistence = alpha + beta + 0.5 * gamma. For EGARCH, persistence is the sum of beta coefficients.

A persistence > 1 implies non-stationarity (IGARCH behavior). Most equity series have persistence in [0.95, 0.999].

See also

garch_fit: Estimate GARCH parameters. garch_forecast: Use persistence for multi-step forecasting.

hawkes_process(events, *, max_time=None, n_iter=500)[source]

Fit a univariate Hawkes (self-exciting) point process via MLE.

The Hawkes process models volatility clustering as a self-exciting point process where each event temporarily increases the intensity (rate) of future events. This is directly analogous to how large price moves tend to cluster in financial markets.

The intensity (conditional rate) is:

\[\lambda(t) = \mu + \sum_{t_i < t} \alpha \exp(-\beta (t - t_i))\]

where \(\mu\) is the background intensity, \(\alpha\) is the jump size per event, and \(\beta\) is the exponential decay rate.

Parameters:

  • events (ndarray) – Array of event times (e.g., times of large returns or trades). Must be sorted in ascending order.

  • max_time (float | None, default: None) – End of observation window. If None, uses max(events) * 1.01.

  • n_iter (int, default: 500) – Maximum iterations for the optimizer. Default 500.

Returns:

  • mu (float) – Background (baseline) intensity.

  • alpha (float) – Excitation magnitude. Larger alpha means each event causes a bigger spike in future event intensity.

  • beta (float) – Decay rate. Larger beta means faster return to baseline intensity after an event.

  • branching_ratio (float) – alpha / beta. Must be < 1 for the process to be stationary. Values near 1 indicate strong self-excitation (volatility clustering).

  • half_life (float) – Time for excitation to decay by 50%: ln(2) / beta.

  • log_likelihood (float) – Maximized log-likelihood.

  • intensity (np.ndarray) – Estimated intensity function evaluated at each event time.

Return type:

dict[str, Any]

Example

>>> import numpy as np
>>> from wraquant.vol.models import hawkes_process
>>> rng = np.random.default_rng(42)
>>> # Simulate clustered events
>>> events = np.sort(rng.exponential(1, 100).cumsum())
>>> result = hawkes_process(events)
>>> result['mu'] > 0
True

Notes

The branching ratio alpha/beta controls the degree of self-excitation. A ratio near 1 means the process is barely stationary and events strongly cluster. In financial markets, branching ratios of 0.5-0.8 are common for trade arrivals.

The Hawkes process is equivalent to an inhomogeneous Poisson process whose intensity depends on the history of past events. It was originally developed for earthquake modeling (Hawkes, 1971) and has been widely applied to high-frequency finance.

Reference: Hawkes, A.G. (1971). “Spectra of some self-exciting and mutually exciting point processes.” Biometrika, 58(1), 83–90.

See also

stochastic_vol_sv: Continuous-time stochastic volatility model. gaussian_mixture_vol: Regime-based volatility clustering.

stochastic_vol_sv(returns, n_particles=500, n_iter=20)[source]

Fit a basic stochastic volatility model via particle filter.

The discrete-time stochastic volatility (SV) model treats log-volatility as a latent AR(1) process:

\[\begin{split}y_t &= \exp(h_t / 2) \, \epsilon_t, \quad \epsilon_t \sim N(0, 1) \\ h_t &= \mu + \phi (h_{t-1} - \mu) + \sigma_\eta \eta_t, \quad \eta_t \sim N(0, 1)\end{split}\]

where \(h_t\) is the log-variance, \(\phi\) controls persistence, \(\mu\) is the long-run mean of log-variance, and \(\sigma_\eta\) is the volatility of volatility.

Parameters:

  • returns (Series | ndarray) – Return series (not prices).

  • n_particles (int, default: 500) – Number of particles for the bootstrap particle filter. More particles = more accurate but slower. Default 500.

  • n_iter (int, default: 20) – Number of iterations for parameter learning via particle marginal MH or iterative filtering. Default 20.

Returns:

  • mu (float) – Long-run mean of log-variance.

  • phi (float) – Persistence of log-variance AR(1). Values near 1 indicate very persistent volatility.

  • sigma_eta (float) – Volatility of volatility.

  • filtered_volatility (np.ndarray) – Filtered estimate of exp(h_t/2), the conditional standard deviation.

  • log_variance (np.ndarray) – Filtered log-variance h_t.

  • log_likelihood (float) – Approximate log-likelihood from the particle filter.

Return type:

dict[str, Any]

Example

>>> import numpy as np
>>> from wraquant.vol.models import stochastic_vol_sv
>>> rng = np.random.default_rng(42)
>>> returns = rng.normal(0, 0.01, 500)
>>> result = stochastic_vol_sv(returns, n_particles=200, n_iter=5)
>>> result['phi'] > 0
True

Notes

The SV model is an alternative to GARCH that treats volatility as an unobserved latent process rather than a deterministic function of past returns. This is arguably more realistic but harder to estimate.

This implementation uses a bootstrap particle filter (sequential Monte Carlo) for likelihood evaluation and a simple grid-based parameter search. For production Bayesian SV estimation, consider PyMC or Stan.

Reference: Taylor, S.J. (1982). “Financial Returns Modelled by the Product of Two Stochastic Processes.” In Time Series Analysis: Theory and Practice 1, O.D. Anderson (ed.), 203–226.

See also

garch_fit: Deterministic volatility filtering (GARCH). gaussian_mixture_vol: Regime-switching volatility. hawkes_process: Point process for event clustering.

gaussian_mixture_vol(returns, n_components=2, random_state=42)[source]

Fit a Gaussian Mixture Model for regime-dependent volatility.

Models the return distribution as a mixture of Gaussian components, each representing a different volatility regime. This captures the empirical observation that returns alternate between calm (low-vol) and turbulent (high-vol) periods.

\[f(r_t) = \sum_{k=1}^{K} \pi_k \, \mathcal{N}(r_t; \mu_k, \sigma_k^2)\]
Parameters:

  • returns (Series | ndarray) – Return series (not prices).

  • n_components (int, default: 2) – Number of Gaussian components (regimes). Default 2 (low-vol and high-vol). Use 3 for an additional crisis regime.

  • random_state (int, default: 42) – Random seed for reproducibility.

Returns:

  • weights (np.ndarray) – Mixing weights (pi_k) for each component, shape (n_components,).

  • means (np.ndarray) – Mean return in each regime.

  • volatilities (np.ndarray) – Volatility (std dev) in each regime, sorted from lowest to highest.

  • regime_probabilities (np.ndarray) – Posterior probability of each regime for each observation, shape (T, n_components).

  • regime_labels (np.ndarray) – Most likely regime for each observation, shape (T,).

  • aic (float) – Akaike Information Criterion.

  • bic (float) – Bayesian Information Criterion.

  • model – Fitted sklearn GaussianMixture object.

Return type:

dict[str, Any]

Example

>>> import numpy as np
>>> from wraquant.vol.models import gaussian_mixture_vol
>>> rng = np.random.default_rng(42)
>>> low_vol = rng.normal(0, 0.005, 300)
>>> high_vol = rng.normal(0, 0.02, 200)
>>> returns = np.concatenate([low_vol, high_vol])
>>> result = gaussian_mixture_vol(returns, n_components=2)
>>> len(result['volatilities']) == 2
True

Notes

The GMM approach is simpler than hidden Markov models but does not model transition dynamics between regimes. For Markov regime-switching, see wraquant.regimes.

The number of components can be selected via AIC/BIC by fitting models with different n_components values.

See also

stochastic_vol_sv: Continuous latent volatility process. garch_fit: GARCH-based volatility modeling.

vol_surface_svi(strikes, forward, total_implied_var, time_to_expiry)[source]

Fit the SVI parameterization to an implied volatility smile.

The Stochastic Volatility Inspired (SVI) model of Gatheral (2004) parameterizes the total implied variance as a function of log-moneyness:

\[w(k) = a + b \left( \rho (k - m) + \sqrt{(k - m)^2 + \sigma^2} \right)\]

where \(k = \ln(K/F)\) is log-moneyness, and the five parameters (a, b, rho, m, sigma) control the level, slope, curvature, translation, and minimum variance of the smile.

Parameters:

  • strikes (ndarray) – Array of strike prices.

  • forward (float) – Forward price of the underlying.

  • total_implied_var (ndarray) – Array of total implied variances (IV^2 * T) corresponding to each strike. Same length as strikes.

  • time_to_expiry (float) – Time to expiry in years.

Returns:

  • params (dict) – Fitted SVI parameters: a (level), b (slope), rho (rotation, -1 to 1), m (translation), sigma (smoothing, > 0).

  • fitted_total_var (np.ndarray) – Fitted total implied variance at each strike.

  • fitted_iv (np.ndarray) – Fitted implied volatility at each strike (sqrt(w / T)).

  • residuals (np.ndarray) – Fitting residuals.

  • rmse (float) – Root mean squared error of the fit.

Return type:

dict[str, Any]

Example

>>> import numpy as np
>>> from wraquant.vol.models import vol_surface_svi
>>> strikes = np.array([90, 95, 100, 105, 110], dtype=float)
>>> forward = 100.0
>>> iv = np.array([0.25, 0.22, 0.20, 0.22, 0.25])
>>> T = 0.25
>>> total_var = iv**2 * T
>>> result = vol_surface_svi(strikes, forward, total_var, T)
>>> result['rmse'] < 0.01
True

Notes

The SVI parameterization is widely used by practitioners for interpolating and extrapolating the volatility smile. It naturally produces realistic smile shapes and can be calibrated to satisfy no-arbitrage constraints (Gatheral and Jacquier, 2014).

Key properties:

  • rho < 0: negative skew (typical for equities).

  • b controls the overall slope of the wings.

  • sigma controls the curvature near the money.

Reference: Gatheral, J. (2004). “A Parsimonious Arbitrage-Free Implied Volatility Parameterization.” Presentation at Global Derivatives & Risk Management, Madrid.

See also

variance_risk_premium: Compute VRP from implied and realized vol.

variance_risk_premium(implied_vol, realized_vol, annualized=True)[source]

Compute the variance risk premium (VRP).

The variance risk premium is the difference between risk-neutral (implied) and physical (realized) variance. A positive VRP indicates that investors pay a premium for volatility protection, which is the empirical norm for equity indices.

\[VRP_t = IV_t^2 - RV_t^2\]

A consistently positive VRP is the economic basis for volatility selling strategies (e.g., selling straddles, variance swaps).

Parameters:

  • implied_vol (Series | ndarray) – Implied volatility series (e.g., VIX / 100 for S&P 500). Should be in the same units as realized_vol (both annualized or both not).

  • realized_vol (Series | ndarray) – Realized volatility series. Same length and alignment as implied_vol.

  • annualized (bool, default: True) – Whether the input volatilities are annualized. If True, the VRP is in annualized variance units.

Returns:

  • vrp (np.ndarray) – Variance risk premium time series. Positive values = implied > realized (normal).

  • mean_vrp (float) – Average VRP over the sample.

  • vrp_ratio (np.ndarray) – IV / RV ratio. Values > 1 indicate implied exceeds realized.

  • pct_positive (float) – Fraction of observations where VRP > 0. Typically 60–80% for equity indices.

  • vol_spread (np.ndarray) – Simple difference IV - RV (in volatility, not variance units).

Return type:

dict[str, Any]

Example

>>> import numpy as np
>>> from wraquant.vol.models import variance_risk_premium
>>> iv = np.array([0.20, 0.22, 0.18, 0.25, 0.19])
>>> rv = np.array([0.15, 0.17, 0.16, 0.20, 0.14])
>>> result = variance_risk_premium(iv, rv)
>>> result['mean_vrp'] > 0
True

Notes

The VRP is one of the most robust risk premiums in financial markets. It reflects the cost of insurance against volatility spikes and is related to jump risk, model uncertainty, and aggregate risk aversion.

Reference: Bollerslev, T., Tauchen, G., and Zhou, H. (2009). “Expected Stock Returns and Variance Risk Premia.” Review of Financial Studies, 22(11), 4463–4492.

See also

garch_fit: Estimate realized conditional volatility. vol_surface_svi: Implied volatility surface fitting.

Realized Volatility

Realized volatility estimators.

Provides various volatility estimators from OHLCV data, including classical, range-based, and high-frequency estimators.

realized_volatility(returns, window=20, annualize=True, periods_per_year=252)[source]

Rolling realized volatility from return series.

Parameters:

  • returns (Series) – Return series.

  • window (int, default: 20) – Rolling window size.

  • annualize (bool, default: True) – Whether to annualize the volatility.

  • periods_per_year (int, default: 252) – Periods per year for annualization.

Return type:

Series

Returns:

Rolling realized volatility series.

parkinson(high, low, window=20, annualize=True, periods_per_year=252)[source]

Parkinson (1980) range-based volatility estimator.

More efficient than close-to-close for continuous processes.

Parameters:

  • high (Series) – High price series.

  • low (Series) – Low price series.

  • window (int, default: 20) – Rolling window size.

  • annualize (bool, default: True) – Whether to annualize.

  • periods_per_year (int, default: 252) – Periods per year.

Return type:

Series

Returns:

Parkinson volatility series.

garman_klass(open_, high, low, close, window=20, annualize=True, periods_per_year=252)[source]

Garman-Klass (1980) volatility estimator.

Uses open, high, low, close for higher efficiency than Parkinson or close-to-close.

Parameters:

  • open – Open price series.

  • high (Series) – High price series.

  • low (Series) – Low price series.

  • close (Series) – Close price series.

  • window (int, default: 20) – Rolling window size.

  • annualize (bool, default: True) – Whether to annualize.

  • periods_per_year (int, default: 252) – Periods per year.

  • open_ (Series)

Return type:

Series

Returns:

Garman-Klass volatility series.

rogers_satchell(open_, high, low, close, window=20, annualize=True, periods_per_year=252)[source]

Rogers-Satchell (1991) volatility estimator.

Handles drift in the price process, unlike Parkinson or Garman-Klass.

Parameters:

  • open – Open price series.

  • high (Series) – High price series.

  • low (Series) – Low price series.

  • close (Series) – Close price series.

  • window (int, default: 20) – Rolling window size.

  • annualize (bool, default: True) – Whether to annualize.

  • periods_per_year (int, default: 252) – Periods per year.

  • open_ (Series)

Return type:

Series

Returns:

Rogers-Satchell volatility series.

yang_zhang(open_, high, low, close, window=20, annualize=True, periods_per_year=252)[source]

Yang-Zhang (2000) volatility estimator.

Combines overnight and Rogers-Satchell estimators for minimum variance under drift and opening jumps.

Parameters:

  • open – Open price series.

  • high (Series) – High price series.

  • low (Series) – Low price series.

  • close (Series) – Close price series.

  • window (int, default: 20) – Rolling window size.

  • annualize (bool, default: True) – Whether to annualize.

  • periods_per_year (int, default: 252) – Periods per year.

  • open_ (Series)

Return type:

Series

Returns:

Yang-Zhang volatility series.

bipower_variation(returns, window=20, annualize=True, periods_per_year=252)[source]

Bipower variation – jump-robust realised volatility estimator.

The bipower variation (BPV) of Barndorff-Nielsen and Shephard (2004) estimates the continuous (diffusive) component of quadratic variation by using the product of adjacent absolute returns instead of squared returns. Because jumps affect only isolated observations, taking the product of consecutive absolute returns washes out the jump contribution asymptotically:

\[BPV_t = \frac{\pi}{2} \sum_{i=2}^{n} |r_{t,i}| \cdot |r_{t,i-1}|\]

where the \(\pi / 2\) factor corrects for the bias introduced by using absolute values of normal random variables (\(E[|z|] = \sqrt{2/\pi}\) for standard normal z).

Parameters:

  • returns (Series) – Return series (not prices).

  • window (int, default: 20) – Rolling window size. Each BPV estimate uses window observations.

  • annualize (bool, default: True) – Whether to annualize by sqrt(periods_per_year).

  • periods_per_year (int, default: 252) – Trading periods per year (252 for daily).

Return type:

Series

Returns:

Rolling bipower variation series (as volatility, i.e. the square root of the bipower variation divided by the window) with the same index as the input.

Example

>>> import pandas as pd, numpy as np
>>> rng = np.random.default_rng(42)
>>> returns = pd.Series(rng.normal(0, 0.01, 200))
>>> bpv = bipower_variation(returns, window=20)
>>> (bpv.dropna() > 0).all()
True

Notes

In the absence of jumps, BPV converges to the same value as standard realised variance. When jumps are present, BPV < RV and the difference RV - BPV estimates the jump component.

This is a rolling window version that computes BPV over a trailing window of window observations at each point.

Reference: Barndorff-Nielsen, O.E. and Shephard, N. (2004). “Power and Bipower Variation with Stochastic Volatility and Jumps.” Journal of Financial Econometrics, 2(1), 1–37.

See also

realized_volatility: Standard rolling realised volatility. jump_test_bns: Formal jump detection test using BPV.

jump_test_bns(returns, window=20, alpha=0.05)[source]

Barndorff-Nielsen & Shephard (2004) jump detection test.

Tests whether jumps are present in a return series by comparing realised variance (RV) to bipower variation (BPV). Under the null hypothesis of no jumps the two estimators are asymptotically equivalent. A significantly positive difference indicates the presence of jumps.

The test statistic is:

\[z = \frac{RV - BPV}{\sqrt{\theta \cdot QP}}\]

where QP is the realised quarticity (estimated via the quad-power quarticity estimator) and theta is a finite-sample correction factor:

\[\theta = \left(\frac{\pi^2}{4} + \pi - 5\right) \cdot \frac{1}{n}\]

The statistic is asymptotically standard normal under the null of no jumps.

Parameters:

  • returns (Series) – Return series (not prices).

  • window (int, default: 20) – Number of observations to use for the test. The test is applied to the most recent window observations. Use the full series length by setting window=len(returns).

  • alpha (float, default: 0.05) – Significance level for the jump detection decision. Default 0.05.

Returns:

  • rv (float) – Realised variance over the window.

  • bpv (float) – Bipower variation over the window.

  • jump_component (float) – max(RV - BPV, 0).

  • continuous_component (float) – BPV.

  • z_statistic (float) – Test statistic (standard normal under the null of no jumps).

  • p_value (float) – One-sided p-value.

  • jump_detected (bool) – True if p-value < alpha.

Return type:

dict[str, object]

Example

>>> import pandas as pd, numpy as np
>>> rng = np.random.default_rng(42)
>>> returns = pd.Series(rng.normal(0, 0.01, 100))
>>> result = jump_test_bns(returns, window=100)
>>> 'z_statistic' in result and 'p_value' in result
True

Notes

The test has low power against small or infrequent jumps in short samples. For daily data, windows of at least 60–100 observations are recommended.

Reference: Barndorff-Nielsen, O.E. and Shephard, N. (2004). “Power and Bipower Variation with Stochastic Volatility and Jumps.” Journal of Financial Econometrics, 2(1), 1–37.

See also

bipower_variation: The jump-robust volatility estimator. realized_volatility: Standard realised volatility.

two_scale_realized_variance(returns, window=20, n_slow=5, annualize=True, periods_per_year=252)[source]

Two-Scale Realised Variance (TSRV) – noise-robust RV estimator.

The TSRV of Zhang, Mykland, and Ait-Sahalia (2005) eliminates the bias caused by market microstructure noise (bid-ask bounce, discrete prices) by combining realised variances computed at two different sampling frequencies. The fast-scale estimator is contaminated by noise; the slow-scale estimator has less noise but more discretisation error. Their optimal linear combination cancels the noise bias:

\[TSRV = RV^{(slow)} - \frac{\bar{n}}{n} RV^{(all)}\]

where \(RV^{(slow)}\) is the average of sub-sampled RVs at a slower frequency, \(RV^{(all)}\) uses all ticks, and \(\bar{n}\) is the average sub-sample size.

Parameters:

  • returns (Series) – Return series. For intraday data, these should be tick-by-tick or high-frequency log returns.

  • window (int, default: 20) – Rolling window size.

  • n_slow (int, default: 5) – Sub-sampling factor for the slow scale. The fast scale uses every observation; the slow scale uses every n_slow-th. Larger values reduce noise but increase discretisation error. Default 5.

  • annualize (bool, default: True) – Whether to annualize.

  • periods_per_year (int, default: 252) – Trading periods per year.

Return type:

Series

Returns:

Rolling TSRV series (as volatility, i.e. square root).

Example

>>> import pandas as pd, numpy as np
>>> rng = np.random.default_rng(42)
>>> returns = pd.Series(rng.normal(0, 0.01, 200))
>>> tsrv = two_scale_realized_variance(returns, window=20)
>>> (tsrv.dropna() > 0).all()
True

Notes

TSRV is consistent for the integrated variance even in the presence of i.i.d. microstructure noise, unlike standard RV which diverges as the sampling frequency increases. It is particularly useful for intraday data sampled at very high frequencies (seconds or ticks).

Reference: Zhang, L., Mykland, P.A., and Ait-Sahalia, Y. (2005). “A Tale of Two Time Scales: Determining Integrated Volatility with Noisy High-Frequency Data.” Journal of the American Statistical Association, 100(472), 1394–1411.

See also

realized_volatility: Standard RV (no noise correction). realized_kernel: Alternative noise-robust estimator.

realized_kernel(returns, window=20, kernel='parzen', annualize=True, periods_per_year=252)[source]

Realised kernel estimator – noise-robust flat-top kernel RV.

The realised kernel of Barndorff-Nielsen, Hansen, Lunde, and Shephard (2008) uses a kernel-weighted sum of autocovariances of returns to produce a consistent, non-negative estimator of integrated variance in the presence of market microstructure noise:

\[RK = \sum_{h=-H}^{H} k\!\left(\frac{h}{H+1}\right) \hat{\gamma}(h)\]

where \(\hat{\gamma}(h) = \sum_i r_i r_{i-h}\) is the sample autocovariance at lag h, and \(k(\cdot)\) is a kernel function with \(k(0) = 1\), \(k(x) = 0\) for \(|x| > 1\).

The Parzen kernel is used by default as it guarantees non-negativity:

\[\begin{split}k(x) = \begin{cases} 1 - 6x^2 + 6|x|^3 & \text{if } |x| \le 1/2 \\ 2(1-|x|)^3 & \text{if } 1/2 < |x| \le 1 \\ 0 & \text{otherwise} \end{cases}\end{split}\]
Parameters:

  • returns (Series) – Return series.

  • window (int, default: 20) – Rolling window size.

  • kernel (str, default: 'parzen') – Kernel function. Currently supported: "parzen" (default). The Parzen kernel guarantees non-negativity.

  • annualize (bool, default: True) – Whether to annualize.

  • periods_per_year (int, default: 252) – Trading periods per year.

Return type:

Series

Returns:

Rolling realised kernel volatility series (square root of the kernel-weighted variance).

Example

>>> import pandas as pd, numpy as np
>>> rng = np.random.default_rng(42)
>>> returns = pd.Series(rng.normal(0, 0.01, 200))
>>> rk = realized_kernel(returns, window=20)
>>> (rk.dropna() > 0).all()
True

Notes

The realised kernel is rate-optimal under i.i.d. noise and achieves a convergence rate of \(n^{-1/5}\) (vs \(n^{-1/4}\) for standard RV with optimal sparse sampling).

The bandwidth H is chosen automatically as \(H = \lceil c \cdot n^{3/5} \rceil\) where c is a constant typically set to 1.

Reference: Barndorff-Nielsen, O.E., Hansen, P.R., Lunde, A., and Shephard, N. (2008). “Designing Realized Kernels to Measure the Ex-Post Variation of Equity Prices in the Presence of Noise.” Econometrica, 76(6), 1481–1536.

See also

two_scale_realized_variance: Alternative noise-robust estimator. realized_volatility: Standard RV without noise correction. bipower_variation: Jump-robust estimator.