Volatility clustering and GARCH models

Take any equity index, record the daily percentage return for several years, and plot the time series. Two things stand out immediately:

1. The returns look roughly centred on zero — there is no obvious upward or downward drift day-to-day.
2. The size of the moves is not constant. Periods of large swings (±2%, ±3%) are followed by more large swings. Periods of calm (±0.2%) are followed by more calm.

This is volatility clustering: the magnitude of today's return is positively correlated with the magnitude of yesterday's return.

It is important to note what clustering does not mean. It does not mean "if the market fell today it will fall again tomorrow." The direction is unpredictable. Only the size clusters.

A simple model treats every daily return as an independent draw from a fixed distribution (say, Gaussian with mean 0 and standard deviation 1%). Under that model, a 3% move and a 0.2% move are equally likely on any given day — volatility is constant.

The real data says something different: if yesterday's move was 3%, today is much more likely to be another big move than a quiet 0.2% day. A model that ignores this will have tails that are far too thin: it will underestimate the probability of back-to-back large moves, which is exactly when risk explodes.

Think of HTTP request rate on a busy service. When you are experiencing a traffic spike, the next 5-minute bucket is more likely to also be high than to have returned to baseline. When traffic is low, it tends to stay low. A naive model that treats each bucket as independent — sampling from the same distribution regardless of recent history — would produce a realistic average request rate but would badly underestimate the probability of sustained spikes, which is exactly what stresses your infrastructure.

Volatility clustering in financial returns is the same phenomenon. The underlying "request generator" (the market) has memory: its current state affects the distribution of near-future outputs.

GARCH (Generalised AutoRegressive Conditional Heteroskedasticity, Engle 1982 / Bollerslev 1986) is the standard model for volatility clustering. The variance on day t is:

sigma^2_t = omega + alpha * epsilon^2_{t-1} + beta * sigma^2_{t-1}

where:

Stationarity (finite long-run variance) requires alpha + beta < 1.

In code, GARCH(1,1) is:

def garch11(returns, omega, alpha, beta, var_init):
    var = var_init
    vars_out = []
    for r in returns:
        var = omega + alpha * r**2 + beta * var
        vars_out.append(var)
    return vars_out

Compare this to an exponential weighted moving average (EWMA) on squared returns (the RiskMetrics model):

def ewma_var(returns, lam, var_init):
    var = var_init
    for r in returns:
        var = (1 - lam) * r**2 + lam * var
    return var

GARCH is EWMA with two differences: (1) an additive floor omega that prevents variance from collapsing to zero even after a long quiet period, and (2) separate alpha and beta weights that do not have to sum to 1 (though they must sum to less than 1 for stationarity).

• Large beta (near 1): variance is persistent. A shock today will still be visible in variance 10 or 20 days later. Real equity markets typically have beta around 0.85–0.95.
• Large alpha: variance reacts sharply to recent squared returns. A single large move immediately inflates tomorrow's variance. Typical equity values are 0.05–0.15.
• alpha + beta near 1: the model is "near-integrated" — shocks are very slow to decay. In the GJR-GARCH paper this is the hard corner of the parameter space (see * Stationarity).
• Small omega: the long-run variance floor is low; all the action is in the alpha and beta terms.

GARCH is symmetric. A +3% return and a -3% return produce exactly the same value of epsilon^2_{t-1} (9%), and therefore exactly the same increase in tomorrow's variance.

Empirically, this is wrong. In equity markets, bad news increases volatility more than equally-sized good news. A -3% crash day raises vol more than a +3% rally day. GARCH cannot capture this.

GJR-GARCH (Glosten, Jagannathan, and Runkle 1993) adds one term that fires only on negative returns:

sigma^2_t = omega + alpha * epsilon^2_{t-1}
                  + gamma * epsilon^2_{t-1} * I{epsilon_{t-1} < 0}
                  + beta * sigma^2_{t-1}

I{epsilon_{t-1} < 0} is an indicator that equals 1 when yesterday's return was negative, and 0 otherwise.

Effect on tomorrow's variance:

When gamma > 0, a negative return amplifies variance more than a positive return of the same size. This is the leverage effect.

The name comes from corporate finance. When a stock price drops, the company's equity-to-debt ratio falls — the same amount of debt is now a larger fraction of the firm's value. Higher financial leverage means higher equity risk, hence higher volatility. In practice, the leverage effect is observed empirically across all equity markets, with or without a clean balance-sheet story behind it.

def gjr_garch_step(prev_return, prev_var, omega, alpha, gamma, beta):
    is_negative = 1.0 if prev_return < 0.0 else 0.0
    asymmetric_alpha = alpha + gamma * is_negative
    return omega + asymmetric_alpha * prev_return**2 + beta * prev_var

That is GJR-GARCH in full. The rest of the model (innovations, simulation loop, parameter estimation) is bookkeeping around this core update.

Under GARCH(1,1), if alpha + beta < 1, the variance process is covariance-stationary: it has a finite long-run (unconditional) average value. That unconditional variance is:

sigma^2_inf = omega / (1 - alpha - beta)

You can think of this as the equilibrium level that variance reverts to after a shock has decayed away.

For GJR-GARCH the stationarity condition is:

kappa = alpha + beta + gamma * p_minus < 1

where p_minus = E[z^2 * I{z < 0}] is the expected squared negative innovation. For a symmetric distribution (Gaussian or symmetric Student-t), p_minus = 0.5. For the skewed distributions used in practice, p_minus differs slightly from 0.5 and depends on the shape parameters.

The unconditional variance under GJR-GARCH is:

sigma^2_inf = omega / (1 - kappa)

As kappa -> 1, the denominator 1 - kappa -> 0 and the unconditional variance diverges. Shocks are so persistent that they never fully die out. This is the "near-integrated" or "near-unit-root" regime.

It is the hardest corner of the parameter space for any numerical method — including the MDN surrogate in the sprint-21 paper — because small changes in parameters produce large changes in the long-run distribution. The paper handles this by engineering a special input feature log(1.01 - kappa) that stretches the near-boundary region and reduces surrogate error there by roughly half.

In a GARCH model, the daily return epsilon_t is modelled as:

epsilon_t = sigma_t * z_t

where z_t is a standardised innovation drawn independently each day.

The simplest choice is z_t ~ N(0, 1) (standard Gaussian). But real equity returns have two systematic departures from Gaussian:

1. Fat tails: extreme moves (more than 3 standard deviations) happen far more often than a Gaussian predicts. A −5% single-day move in a broad equity index is a "5-sigma event" under Gaussian assumptions, implying it should happen once in millions of years. In reality it happens every few years.
2. Negative skew: large negative moves are more extreme than large positive moves. Equity distributions have a longer left tail.

The GJR-GARCH NN paper uses a skewed Student-t distribution for z_t, parameterised by:

For large nu (say, nu > 30), the skewed t converges to a Gaussian with skewness lambda. For nu = 4, the tails are substantially fatter than Gaussian — a good match for daily equity data.

Because the skewed t is not symmetric, p_minus = E[z^2 * I{z < 0}] is no longer exactly 0.5. It depends on both nu and lambda and must be computed numerically from the distribution's CDF. This feeds directly into the persistence coefficient kappa and the unconditional variance formula.

The innovation z_t is not a Gaussian or a skewed t draw. It is a draw from the feed's own Gaussian mixture, the GMM of The GMM component, standardised to mean zero and variance one:

z_t = (x_t - m) / s       x_t drawn from the mixture, m its mean, s its standard deviation
epsilon_t = sigma_t * z_t

The two models then do separate work. The mixture gives the shape of a tick's shock: its fat tails, its asymmetry and the regimes its components stand for. The recursion gives the size of the shock as a function of the days before, which an independent mixture draw cannot. The skewed Student-t of Hansen (1994) would give the shape with two parameters; the mixture is used because its components are already part of each FX spot configuration and are already calibrated from public fixings.

For a mixture, p_minus = E[z^2 * I{z < 0}] has a closed form. Each component k, with weight w_k, mean mu_k and standard deviation s_k after standardisation, contributes the second moment of a normal truncated at zero:

p_minus = sum_k w_k * [ (mu_k^2 + s_k^2) * Phi(-mu_k / s_k) - mu_k * s_k * phi(mu_k / s_k) ]

where Phi and phi are the standard normal distribution and density.

The recursion reacts to a day's shock, which is the sum of the day's n ticks and not one tick. Summing independent draws keeps the variance additive and thins the skewness s to s / sqrt(n), so a day's shock is closer to normal than a tick's. The first Edgeworth term of the sum gives its p_minus:

p_minus_day = 1/2 - s / (3 * sqrt(2 * pi) * sqrt(n))

which is one half for a symmetric mixture and tends to one half as n grows. A day of one tick takes the closed form above. The persistence kappa and the unconditional variance follow from p_minus_day as in * Stationarity. Taking the tick's p_minus instead overstates the asymmetry's share of the persistence for a skewed mixture, so variance targeting sets omega too low and the daily variance settles below its intended level.

In equities a fall is bad news, so gamma > 0. In FX a fall of the pair is a rise of its base currency, and which direction raises volatility depends on the pair: a pair whose quote currency is a safe haven tends to fall in a crisis, a pair whose base currency is one tends to rise. So gamma may take either sign. The variance stays positive if

omega > 0,   alpha >= 0,   beta >= 0,   alpha + gamma >= 0

and a negative gamma, bounded by -alpha, makes a rise of the pair raise variance more than a fall of the same size.

The parameters are estimated from the daily fixings of the pair, five years of which every environment holds; see Synthetic prices are calibrated from free public sources.

1. Variance targeting. The unconditional variance is set to the sample variance of the daily log returns, and omega = sigma^2_inf * (1 - kappa). One parameter leaves the fit, and the long-run level is tied to the data directly (Engle and Mezrich, 1996).
2. Quasi maximum likelihood. alpha, gamma and beta maximise the Gaussian likelihood of the returns under the bounds above. The likelihood is Gaussian although the innovation is a mixture: the estimator is consistent when the conditional mean and variance are correctly specified, whatever the shape of the innovation (Bollerslev and Wooldridge, 1992).
3. The starting variance. sigma^2_0 is the last value of the fitted variance filtered through the whole history, so a feed that starts in a turbulent period starts turbulent.
4. Authored values. A pair with too little history for a stable fit takes an authored shape, as it takes authored mixture components. An authored shape is a calibration like a fitted one, so it is reviewed and recorded the same way.

A feed never fits its own variance. The estimate above is a calibration with three stages, each stored as a versioned record, so every number a feed runs traces back to the fixings it came from and to the person who accepted it. The first two stages run inside the class calibration run, beside the level and the spread, and the third is a review; Calibrating a synthetic feed: from fixings to a reviewed run draws the whole flow.

1. Capture the inputs. A capture reads the pair's fixings under the asset class's calibration policy, which states the variance window, the order of the providers and whether an inverse may be read. It stores a variance calibration input: the pair, the item and provider read, the policy and its version, the window, whether the item was inverted, and one point per fixing date. A point holds the fixing as stored, the return it gives in the pair's market convention, and whether it was used; a point left out, such as a value that is not a number, keeps its reason. Nothing is fitted at this stage.
2. Fit and review. A fit reads one stored input, never the live fixings, so the same input always gives the same result. It stores a variance calibration that names its input and the method's version, and holds what a reviewer needs to judge it: alpha, gamma, beta, the persistence, the long-run and starting variance, the likelihood, and whether the starts of the simplex agreed. A calibration is proposed when it is stored. A person approves or rejects it and states why; approving one supersedes the pair's earlier approval, so one approved calibration is current. A fit that cannot be made, for too few returns or a shape at the stationarity bound, is stored as rejected with its reason, so the attempt is kept.
3. Apply. A feed's configuration names the approved calibration it runs. Choosing another is a new version of the configuration, so the feed's history shows when each calibration came into use. When a run starts it reads the calibration it names and fits nothing, and its start record names that calibration and its version. A feed whose calibration is no longer approved does not start, and says why.

The chain reads in both directions: from a run to its calibration, to the input the calibration was fitted from, to the fixings of that input; and from a calibration to the feeds that run it. A negative gamma, which says that falls calm the market rather than excite it, is the kind of estimate this review exists to catch before a feed runs it.

The recursion is daily, and a feed ticks at gaps its arrival process draws, as Tick arrival describes. The market day ends at the FX day roll, 17:00 New York time, on the market's operational clock, and sigma^2_t holds for the whole day. A tick that covers a fraction of the day takes that fraction of sigma^2_t as its variance: the mixture's deviations from its mean are scaled so the mixture's daily variance is sigma^2_t, and its drift is left alone. At the roll, epsilon_t is the log return of the generated path over the day. Saturday and Sunday are not market days, so the variance does not change across them.

The model is native to ores.analytics.quant, with no QuantLib code. math::gjr_garch_model, in ores.analytics.quant/math/gjr_garch.hpp, can only be built by make(), which refuses a shape outside the admissible region, so a model in hand always has a positive variance and a finite long-run level; omega follows from the long-run variance. The same header holds the filter, the Gaussian cost, the mixture's closed-form p_minus (negative_second_moment), its daily counterpart (daily_negative_second_moment) and fit_gjr_garch, which fits from several starts by the native simplex of math/nelder_mead.hpp and keeps the best. service::gjr_garch_variance_process wraps any process that honours set_volatility_multiplier(), which the two mixture processes do, and process_factory::make_gjr_garch_process builds it with the mixture's own daily variance as the long-run level. The decorator reads the day's return once, at the roll, so it adds a few nanoseconds a tick to the mixture's step. The FX spot feed builds the process from the approved calibration its configuration names, through calibrated_variance in ores.synthetic.api/feeds/variance_calibration.hpp, which refuses a calibration that is missing, of another pair, not approved or shapeless. The same header holds capture_variance_points and fit_variance, which the class calibration run, synthetic.v1.ops.calibrate_feeds, calls on the policy's variance window of each GJR-GARCH feed's item; fit_variance calls fit_gjr_garch on the points of one stored input, and the fit reports the starts it tried and whether they agreed. The review is synthetic.v1.ops.review_variance_calibration.

The tests reproduce QuantLib's GARCH calculation test values, compare the cost with QuantLib's Garch11 cost function ported as an independent expectation, and repeat QuantLib's calibration and initial-guess tests with this project's random numbers: the fit recovers the parameters that drew a series, and every start reaches the same optimum.

Volatility models split into two broad families by time structure:

• Discrete-time models (GARCH family): variance is updated once per day using a difference equation. Parameters are estimated from daily return data. Simulation steps are one trading day.
• Continuous-time models (Heston, SABR, rough Bergomi, …): variance follows a stochastic differential equation (SDE) that can be sampled at any frequency. Parameters are calibrated to option prices.

The Heston model (1993) is the continuous-time model most directly analogous to GARCH. In Heston:

• The asset price follows Geometric Brownian Motion.
• The variance v_t follows its own SDE — a mean-reverting CIR process:

dv_t = kappa * (v_bar - v_t) * dt + xi * sqrt(v_t) * dW_v

• The two Brownian motions (price and variance) are correlated with correlation rho.

Heston parameters map onto GARCH concepts as follows:

The key insight: Heston's rho < 0 and GJR-GARCH's gamma > 0 both capture the same empirical phenomenon (negative shocks increase variance more than positive shocks), just in different mathematical frameworks.

QuantLib mainline does not include a GARCH pricing engine. There is a GARCH11 class in the QuantLib-Ext component (an unofficial extension library), which fits GARCH(1,1) parameters to a historical return series and computes the unconditional variance. It is a calibration utility, not a pricing engine.

For option pricing under GARCH dynamics, the standard path is either:

1. Monte Carlo simulation (expensive, noisy).
2. A surrogate model trained on MC output (what the sprint-21 paper does).

For the continuous-time leverage-effect story, QuantLib is fully equipped. HestonModel and HestonProcess are mainline classes. ORE uses Heston natively for equity and FX vol surface calibration. Setting rho to a large negative value (say, -0.7) in a Heston calibration will produce an implied vol surface with a strong negative skew — the same qualitative behaviour as a GJR-GARCH model with large gamma.

The practical takeaway: if you see HestonModel in ORE configuration, it is capturing the leverage effect via rho. If you see GJR-GARCH in sprint-21 documents, it is capturing the same effect via gamma, but in daily discrete time.

A GARCH or GJR-GARCH model does not produce a QuantLib option price directly. The output of a GJR-GARCH simulation is a distribution of terminal returns (or, via the MDN surrogate, a Gaussian mixture approximation of that distribution). To use this as QuantLib input you convert it to an implied vol surface: evaluate the GJR-GARCH option price at each (strike, maturity) point, then back out the Black-Scholes implied volatility via bisection. The resulting grid can be loaded into ORE as a standard EQUITY_OPTION/RATE vol surface.

Van den Berg (arXiv:2606.15502, June 2026) trains a Mixture Density Network (MDN) to replace Monte Carlo simulation of GJR-GARCH paths for option pricing.

The problem is that GJR-GARCH has no closed-form option price. Pricing a single European call under GJR-GARCH dynamics normally requires running millions of simulated paths and averaging the discounted payoff — which takes hundreds of seconds for a full vol surface. The MDN learns to predict the terminal return distribution directly from the GJR-GARCH parameters and maturity, replacing the MC simulation with a network forward pass that takes microseconds.

The paper uses GJR-GARCH in the Q-measure (risk-neutral) setting with a risk premium adjustment. This means the model is calibrated to produce risk-neutral option prices (the prices you would observe in the market), not physical-measure scenario paths. See Probability measures: P (real-world) and Q (risk-neutral) for the P/Q distinction.

The MDN takes seven scalar inputs (in "dimensionless reduced form"):

(alpha, gamma * p_minus, beta, sigma_0'^2, nu, lambda, T)

where T is the option maturity in trading days. These seven numbers fully specify the GJR-GARCH risk-neutral pricing problem. The output is a 128-component Gaussian mixture that approximates the terminal log-return distribution. European option prices follow from a weighted sum of Black-Scholes formulas over the mixture components — a closed-form calculation once the mixture parameters are known.

The MDN is approximately 400,000 times faster than matched-accuracy MC.

Parameters near the stationarity boundary (kappa -> 1, where variance shocks become highly persistent) are the hardest for the MDN. The 99th percentile CDF error is 4–6 times higher near kappa = 1 than in the centre of the parameter space. The paper mitigates this with an engineered input feature log(1.01 - kappa).

ORE Studio's synthetic pipeline does not use the MDN. GJR-GARCH appears in it under the real-world measure, as the daily variance of the synthetic FX spot described in * Daily GJR-GARCH variance of synthetic FX spot. The synthetic volatility surfaces are not computed from the GARCH state: they evolve through SSVI parameter processes of their own, whose shocks are correlated with the spot shocks, as Synthetic market data generation: approach describes. The paper remains the reference for pricing options under GJR-GARCH dynamics, which the pipeline does not need.

• Engle, R. F. (1982). "Autoregressive conditional heteroscedasticity with estimates of the variance of United Kingdom inflation." Econometrica, 50(4), 987-1007. Listed in Engle's bibliography.
• Bollerslev, T. (1986). "Generalized autoregressive conditional heteroskedasticity." Journal of Econometrics, 31(3), 307-327. doi:10.1016/0304-4076(86)90063-1.
• 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. doi:10.1111/j.1540-6261.1993.tb05128.x.
• Bollerslev, T. and Wooldridge, J. M. (1992). "Quasi-maximum likelihood estimation and inference in dynamic models with time-varying covariances." Econometric Reviews, 11(2), 143-172. doi:10.1080/07474939208800229.
• Engle, R. F. and Mezrich, J. (1996). "GARCH for groups." Risk, 9(8), 36-40. EconBiz record. Variance targeting; the article itself not read.
• Hansen, B. E. (1994). "Autoregressive conditional density estimation." International Economic Review, 35(3), 705-730. doi:10.2307/2527081. The skewed Student-t.
• 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. doi:10.1002/jae.1234.