Volatility
Models whose variance is itself a process — GARCH, EGARCH, Heston, regime switching and AR(1) — with their stationarity tests.
Prelude.Volatility covers the models in which the variance of a series changes over time: GARCH(1, 1) and EGARCH for returns, the Heston stochastic-volatility model, Markov regime switching, and the first-order autoregression they are usually compared with. Each model is given as a one-step update you iterate yourself, plus the closed-form quantities that decide how it behaves.
The phenomenon these models exist to reproduce is volatility clustering: returns are nearly uncorrelated, but their squares are strongly correlated — large moves follow large moves, whatever their sign. A model can have the right average variance and still get this completely wrong, so it is the property worth checking in a fitted model.
For GARCH the number to read first is the persistence \(\alpha + \beta\). Below one, the variance has a long-run level to return to and shocks decay geometrically; at one there is no such level and every shock is permanent. A model fitted at 0.999 is not “nearly stationary” — it behaves as a different object.
Randomness is never hidden: every update that needs noise takes the draw (shock, z, dw) as an argument, so you generate it in C# from whatever seeded stream you use. Tests such as garchStationary return 1 for true and 0 for false, written below with the Iverson bracket \([\,\cdot\,]\). Divisors are guarded with \(\operatorname{nonzero}(d) = d + 10^{-300}\,[\,\lvert d \rvert < 10^{-300}\,]\) (see Chains), which equals \(d\) everywhere except within \(10^{-300}\) of zero.
open Prelude.Volatility
GARCH(1, 1)
Tomorrow’s variance is a weighted combination of a constant, today’s squared shock and today’s variance. The constant \(\omega\) sets the level, \(\alpha\) says how much a new shock matters, and \(\beta\) how much of the old variance carries over. \(\omega\), \(\alpha\) and \(\beta\) must be non-negative for the variance to stay positive.
\[\sigma^2_{t+1} = \omega + \alpha\,\varepsilon_t^2 + \beta\,\sigma_t^2\]garchNext
garchNext(variance, shock, omega, alpha, beta) · variance today’s variance scalar · shock today’s shock (return less its mean) scalar · omega level scalar · alpha shock weight scalar · beta variance weight scalar → tomorrow’s variance scalar
Let \(\sigma_t^2 \ge 0\) (the argument variance) be today’s variance, \(\varepsilon_t \in \mathbb{R}\) (the argument shock) today’s shock, and let \(\omega, \alpha, \beta \ge 0\) be the level, shock weight and variance weight. Then
One step of the GARCH(1, 1) variance recursion. Iterate it along a return series to filter the conditional variance.
Intuition Tomorrow’s variance mixes a floor \(\omega\), today’s surprise squared (a big move of either sign raises risk) and today’s variance carried forward. After a shock, the expected variance decays back by a factor \(\alpha + \beta\) per step. That is how calm and stormy stretches cluster.
Use it when you filter the conditional variance along a return series for volatility forecasts, value-at-risk, position sizing or option inputs, or inside a likelihood to fit \(\omega\), \(\alpha\) and \(\beta\). Feed it demeaned returns, and start it at the sample variance or at unconditional. If falls should raise volatility more than rises of the same size, use egarchNext.
volatility-forecastingrisk-measurementlikelihood-evaluation
def garchNext(variance, shock, omega, alpha, beta) ≔ omega + alpha · shock · shock + beta · variance
persistence
persistence(alpha, beta) · alpha shock weight scalar · beta variance weight scalar → persistence scalar
Let \(\alpha, \beta \ge 0\) be the shock and variance weights. Then
\[\operatorname{persistence}(\alpha, \beta) = \alpha + \beta .\]The persistence decides whether the model is stationary and how slowly shocks die away.
Intuition The share of today’s excess variance that survives to tomorrow on average. It is the decay factor of a volatility shock: at 0.9 a spike fades within days, and at 0.99 it lingers for months of daily data.
Use it when you read a fitted GARCH model, before anything else. Values from 0.95 to 0.99 are common for daily financial returns, and 1 or more means there is no long-run level. Convert it to periods with shockHalfLife, or test it with garchStationary. Persistence very close to 1 can also be a sign of unmodelled regime shifts, so compare a regime model (regimeDuration).
volatility-forecastingmodel-interpretation
def persistence(alpha, beta) ≔ alpha + beta
garchStationary
garchStationary(alpha, beta) · alpha shock weight scalar · beta variance weight scalar → 1 or 0 scalar
Let \(\alpha, \beta \ge 0\) be the shock and variance weights. Then
\[\operatorname{garchStationary}(\alpha, \beta) = [\,\alpha + \beta < 1\,] .\]It is 1 when the variance has an unconditional level, 0 otherwise.
Intuition The expected variance follows \(E[\sigma^2_{t+1}] = \omega + (\alpha+\beta)\,E[\sigma^2_t]\). With \(\alpha + \beta < 1\) that recursion settles at a fixed level. At 1 or more it grows without bound or keeps every shock forever.
Use it when you validate fitted or proposed parameters before computing unconditional, shockHalfLife or long-horizon forecasts, all of which are meaningless otherwise. It also works as a constraint during estimation. A persistence of 0.999 passes, but it behaves almost like the non-stationary case, so read shockHalfLife as well.
stability-checkvolatility-forecasting
def garchStationary(alpha, beta) ≔ persistence(alpha, beta) < 1
unconditional
unconditional(omega, alpha, beta) · omega level scalar · alpha shock weight scalar · beta variance weight scalar → long-run variance scalar
Let \(\omega \ge 0\) be the level and \(\alpha, \beta \ge 0\) the weights, with \(\alpha + \beta < 1\). Then
\[\operatorname{unconditional}(\omega, \alpha, \beta) = \bar{\sigma}^2 = \frac{\omega}{\operatorname{nonzero}(1 - \alpha - \beta)} .\]The unconditional variance the process returns to. Meaningful only for a stationary model.
Intuition The fixed point of the expected-variance recursion: the level at which the \(\omega\) added each step exactly replaces the share \(1 - \alpha - \beta\) that decays. Forecasts far ahead converge to it.
Use it when you need long-horizon volatility (annual risk, inputs for long-dated options), a starting value for the filter, or a check on a fit, since it should be close to the sample variance of the returns. For variance targeting, fix it at the sample variance, set \(\omega\) to that variance times \(1 - \alpha - \beta\), and fit only \(\alpha\) and \(\beta\). It is only meaningful when garchStationary is 1.
volatility-forecastingrisk-measurementparameter-estimation
def unconditional(omega, alpha, beta) ≔ omega / nonzero(1 − persistence(alpha, beta))
shockHalfLife
shockHalfLife(alpha, beta) · alpha shock weight scalar · beta variance weight scalar → half-life in periods scalar
Let \(\alpha, \beta \ge 0\) be the weights, with \(0 < \alpha + \beta < 1\). Then
\[\operatorname{shockHalfLife}(\alpha, \beta) = \frac{\ln 2}{\operatorname{nonzero}\big(-\ln \operatorname{nonzero}(\alpha + \beta)\big)} .\]The number of periods for half of a variance shock to decay. This is what persistence means in practical units: at 0.95 a shock is half gone in about 13.5 periods, at 0.99 in about sixty-nine.
Intuition Persistence in units of time. The excess variance after a shock shrinks by a factor \(\alpha + \beta\) each period, so the half-life is the \(n\) at which \((\alpha+\beta)^n = \tfrac12\).
Use it when you explain or sanity-check a fit: how many trading days a volatility spike lasts, whether risk limits should adapt over days or over months, how two models compare. It is in the periods of your data. It needs \(0 < \alpha + \beta < 1\). At exactly 1 the guard returns an astronomically large number, and above 1 the result is negative. The same idea for other models is halfLife in Diffusions and arHalfLife here.
volatility-forecastingmodel-interpretation
def shockHalfLife(alpha, beta) ≔ ln(2) / nonzero(0 − ln(nonzero(persistence(alpha, beta))))
garchKurtosis
garchKurtosis(alpha, beta) · alpha shock weight scalar · beta variance weight scalar → kurtosis scalar
Let \(\alpha, \beta \ge 0\) be the weights, and write \(p = \alpha + \beta\), with \(1 - p^2 - 2\alpha^2 > 0\). Then
\[\operatorname{garchKurtosis}(\alpha, \beta) = \kappa = \frac{3\,(1 - p^2)}{\operatorname{nonzero}(1 - p^2 - 2\alpha^2)} .\]The unconditional kurtosis of GARCH(1, 1) returns driven by normal innovations. It exceeds 3 whenever \(\alpha > 0\): a varying variance produces fat tails without assuming them. The fourth moment exists only while the denominator is positive; outside that region the number returned has no meaning.
Intuition Even with normal shocks, mixing calm stretches (small variance) with stormy ones (large variance) produces more tiny and more huge returns than a single normal distribution would. Those are fat tails. The larger \(\alpha\), the more the variance itself moves and the fatter the tails.
Use it when you check whether a fitted GARCH explains the kurtosis you measure in returns. If the data are more fat-tailed than this, use fat-tailed innovations or add jumps (compare mertonKurtosis in Diffusions). The fourth moment exists only when \(1 - (\alpha+\beta)^2 - 2\alpha^2 > 0\), and nothing checks that, so test the condition yourself.
risk-measurementmodel-validation
def garchKurtosis(alpha, beta) ≔ 3 · (1 − persistence(alpha, beta) · persistence(alpha, beta)) / nonzero(1 − persistence(alpha, beta) · persistence(alpha, beta) − 2 · alpha · alpha)
EGARCH
EGARCH models the logarithm of the variance, which is free to take any value, so it needs no positivity constraints on its parameters. It also captures the leverage effect — a fall raises volatility more than a rise of the same size — with one extra term \(\gamma\); \(\gamma = 0\) is the symmetric model.
\[\ln \sigma^2_{t+1} = \omega + \beta \ln \sigma^2_t + \alpha\left(|z_t| - \sqrt{2/\pi}\right) + \gamma z_t, \qquad z_t = \varepsilon_t / \sigma_t\]egarchNext
egarchNext(logVariance, shock, omega, alpha, gamma, beta) · logVariance today’s log variance scalar · shock today’s shock scalar · omega level scalar · alpha size effect scalar · gamma sign (leverage) effect scalar · beta persistence scalar → tomorrow’s log variance scalar
Let \(h_t = \ln \sigma_t^2 \in \mathbb{R}\) (the argument logVariance) be today’s log variance, \(\varepsilon_t \in \mathbb{R}\) (the argument shock) today’s shock, and \(\omega, \alpha, \gamma, \beta \in \mathbb{R}\) the parameters. Let \(z_t = \operatorname{standardised}(\varepsilon_t, e^{h_t}) = \varepsilon_t / \sqrt{\max(e^{h_t}, 10^{-300})}\). Then
where \(0.7978845608028654 = \sqrt{2/\pi}\) is the mean of \(\lvert z \rvert\) for a standard normal \(z\). One step of the EGARCH(1, 1) recursion; the shock is standardised by today’s volatility inside the formula.
Intuition Working on the log of the variance keeps the variance positive for any parameter values. The standardised shock enters twice. Its size, measured against the normal average \(\sqrt{2/\pi}\), raises variance after big moves of either sign. Its sign, through \(\gamma z\), lets falls and rises act differently.
Use it when you model equity or index volatility, where falls raise risk more than rallies, or when fitting GARCH keeps running into its non-negativity limits. Iterate on the log variance and exponentiate to get the variance. The centring assumes normal shocks. With \(\gamma = 0\) the model is symmetric. When you need neither asymmetry nor the log form, garchNext is simpler.
volatility-forecastingdownside-riskrisk-measurement
def egarchNext(logVariance, shock, omega, alpha, gamma, beta) ≔ omega + beta · logVariance + alpha · (abs(standardised(shock, exp(logVariance))) − 0.7978845608028654) + gamma · standardised(shock, exp(logVariance))
asymmetric
asymmetric(gamma) · gamma sign effect scalar → 1 or 0 scalar
Let \(\gamma \in \mathbb{R}\) be the sign effect. Then
\[\operatorname{asymmetric}(\gamma) = [\,\gamma < 0\,] .\]It is 1 in the usual leverage direction, in which negative shocks raise volatility more.
Intuition With \(\gamma < 0\), a negative standardised shock adds \(\lvert\gamma z\rvert\) to the log variance and a positive one subtracts the same amount. Falls stoke volatility and rallies calm it.
Use it when you interpret a fitted EGARCH: to confirm the leverage effect typical of equities, or to spot the opposite sign (\(\gamma > 0\)), which some commodity markets show when price spikes bring turbulence. A sign alone can be noise, so check that \(\hat\gamma\) is significantly different from zero before relying on the flag.
downside-riskmodel-interpretation
def asymmetric(gamma) ≔ gamma < 0
egarchStationary
egarchStationary(beta) · beta persistence scalar → 1 or 0 scalar
Let \(\beta \in \mathbb{R}\) be the persistence. Then
\[\operatorname{egarchStationary}(\beta) = [\,\lvert \beta \rvert < 1\,] ,\]the stationarity condition for EGARCH(1, 1).
Intuition In EGARCH the log variance carries into the next step with factor \(\beta\), so, as in an AR(1), it returns to a level only if \(\lvert\beta\rvert < 1\). \(\alpha\) and \(\gamma\) change how shocks enter, not whether they fade.
Use it when you validate EGARCH parameters before long-horizon forecasts. Unlike garchStationary, the test involves \(\beta\) alone. The half-life of a log-variance shock is arHalfLife(beta).
stability-checkvolatility-forecasting
def egarchStationary(beta) ≔ abs(beta) < 1
Building blocks
standardised(shock, variance)— Let \(\varepsilon \in \mathbb{R}\) and \(v \in \mathbb{R}\). Then \(\operatorname{standardised}(\varepsilon, v) = \varepsilon / \operatorname{nonzero}\big(\sqrt{\max(v, 10^{-300})}\big) = \varepsilon / \sqrt{\max(v, 10^{-300})}\), the guard never acting because the root is at least \(10^{-150}\). The shock divided by the volatility, the square root of the variance floored at a tiny positive number.
Heston
In the Heston model the price and its variance each have their own Brownian driver, and the correlation \(\rho\) between them is what makes the model: it is negative in equity markets and produces the skewed volatility smile. Independent drivers would give a symmetric smile and no leverage effect.
\[dv = \kappa(\theta - v)\,dt + \xi\sqrt{v}\,dW^v, \qquad dS = S\left(\mu\,dt + \sqrt{v}\,dW^S\right)\]The steps below are Euler steps with the variance truncated at zero. The increments dw are the caller’s: draw two independent standard normals, correlate them with tied, and scale each by \(\sqrt{dt}\).
tied
tied(first, second, rho) · first a standard normal draw scalar or tensor · second an independent standard normal draw same shape · rho correlation scalar → a draw correlated with first same shape
Let \(z_1, z_2\) (the arguments first and second) be independent standard normal draws of the same shape, taken entrywise, and let \(\rho \in [-1, 1]\) be the correlation. Then
A second normal draw with correlation \(\rho\) to the first — the two-dimensional Cholesky factor, exact.
Intuition Keep a share \(\rho\) of the first draw and fill the rest with fresh independent noise, weighted so the total variance stays 1. The result has correlation exactly \(\rho\) with the first draw.
Use it when you correlate the price and variance drivers of Heston (usually \(\rho < 0\) for equities), the shocks of two assets, or any pair of normals, then scale each with increment. For three or more correlated variables you need a full Cholesky factor (cholesky in Linear algebra). Nothing rejects a \(\rho\) outside \([-1, 1]\), and the result then has variance \(\rho^2 > 1\).
simulationdependence-modellingoption-pricing
def tied(first, second, rho) ≔ rho · first + sqrt(max(1 − rho · rho, 0)) · second
hestonVariance
hestonVariance(v, kappa, theta, xi, dt, dw) · v current variance scalar · kappa mean-reversion speed scalar · theta long-run variance scalar · xi volatility of variance scalar · dt time step scalar · dw Brownian increment for the variance scalar → next variance scalar
Let \(v \in \mathbb{R}\) be the current variance, \(\kappa > 0\) the mean-reversion speed, \(\theta > 0\) the long-run variance, \(\xi > 0\) the volatility of variance, \(\Delta t > 0\) (the argument dt) the time step and \(\Delta W \in \mathbb{R}\) (the argument dw) the Brownian increment. Then
One Euler step of the variance, truncated at zero.
Intuition The variance is itself a square-root diffusion: pulled toward \(\theta\) at speed \(\kappa\), with noise \(\xi\sqrt{v}\) that fades near zero. Clamping each step at zero stops a discrete overshoot from producing a negative variance the price step cannot use.
Use it when you simulate stochastic volatility to price options with a smile or skew, build scenarios with volatility spikes, or test hedging. Give it a dw correlated with the price increment through tied. Euler with a clamp is biased for large steps, especially when hestonFeller fails, so keep dt small. The clamp absorbs the state at zero, which is slightly different from fullTruncation in Diffusions, where the unclamped state is kept.
option-pricingsimulationnumerical-stability
def hestonVariance(v, kappa, theta, xi, dt, dw) ≔ max(v + kappa · (theta − v) · dt + xi · sqrt(max(v, 0)) · dw, 0)
hestonPrice
hestonPrice(x, v, mu, dt, dw) · x current price scalar · v current variance scalar · mu drift scalar · dt time step scalar · dw Brownian increment for the price scalar → next price scalar
Let \(x > 0\) be the current price, \(v \in \mathbb{R}\) the current variance, \(\mu \in \mathbb{R}\) the drift, \(\Delta t > 0\) (the argument dt) the time step and \(\Delta W \in \mathbb{R}\) (the argument dw) the Brownian increment for the price. Then
One Euler step of the price.
Intuition A geometric Brownian step whose volatility is read off the current variance path: \(\sqrt{v}\) takes the place of a constant \(\sigma\). Calm stretches move the price a little and stormy ones a lot.
Use it when you advance the price alongside hestonVariance in a Monte Carlo run, to price path-dependent or out-of-the-money options or to measure risk under changing volatility. Use the variance from the start of the step. It is an Euler step on the price, so a very large \(\sqrt{v}\,\Delta W\) can make the price negative. A log-price step, \(x\,e^{(\mu - v/2)\Delta t + \sqrt{v}\,\Delta W}\), avoids that if you write it yourself.
option-pricingasset-pricingsimulation
def hestonPrice(x, v, mu, dt, dw) ≔ x · (1 + mu · dt + sqrt(max(v, 0)) · dw)
hestonFeller
hestonFeller(kappa, theta, xi) · kappa mean-reversion speed scalar · theta long-run variance scalar · xi volatility of variance scalar → 1 or 0 scalar
Let \(\kappa, \theta, \xi > 0\) be the mean-reversion speed, long-run variance and volatility of variance. Then
\[\operatorname{hestonFeller}(\kappa, \theta, \xi) = [\,2\kappa\theta \ge \xi^2\,] .\]It is 1 when the Feller condition holds, under which the continuous-time variance never reaches zero.
Intuition Near zero, the upward pull on the variance (about \(\kappa\theta\)) must be strong enough relative to the volatility of variance (\(\xi^2\)) for the variance to stay strictly positive.
Use it when you check calibrated Heston parameters before simulating. A 0 does not make the model wrong, and parameters calibrated to equity smiles often fail the test. It does mean variance paths touch zero often, so the bias from clamping grows and dt should shrink. It is fellerHolds in Diffusions, with \(\kappa, \theta, \xi\) in place of \(\theta, \mu, \sigma\).
stability-checknumerical-stabilityoption-pricing
def hestonFeller(kappa, theta, xi) ≔ 2 · kappa · theta ≥ xi · xi
Markov regime switching
A regime-switching model is a mixture whose weights follow a Markov chain, so regimes persist: a calm market stays calm for a while and a turbulent one stays turbulent. transitions is the row-stochastic regime transition matrix, and at is a probability vector over regimes (for example the current filtered regime probabilities).
regimeDuration
regimeDuration(transitions) · transitions regime transition matrix matrix s×s → expected duration of each regime vector s
Let \(P \in [0, 1]^{s \times s}\) (the argument transitions) be a row-stochastic transition matrix. Then, for \(i = 1, \dots, s\),
The expected number of periods spent in each regime once entered.
Intuition Each period the regime stays with probability \(p_{ii}\), so the length of a stay is geometric with mean \(1/(1 - p_{ii})\): 0.9 means about ten periods, 0.99 about a hundred.
Use it when you interpret a fitted regime-switching model, for example how long bull and bear markets, recessions or high-volatility spells last, and check that the estimates are plausible. The durations are in the period of the data the matrix was fitted on. An absorbing regime (\(p_{ii} = 1\)) returns an astronomically large number, not infinity.
regime-modellingmodel-interpretation
def regimeDuration(transitions) ≔ 1 / nonzero(1 − diag(transitions))
regimeMean
regimeMean(at, means) · at regime probabilities vector s · means mean in each regime vector s → mixture mean scalar
Let \(p \in \Delta^{s-1}\) (the argument at) be the regime probabilities and \(\mu \in \mathbb{R}^s\) (the argument means) the regime means. Then
The probability-weighted mean.
Intuition When you are unsure which regime you are in, the expected value is each regime’s mean weighted by how likely that regime is.
Use it when you forecast next-period returns or growth from filtered regime probabilities, or report a blended expectation. To forecast \(k\) periods ahead, first push at through the transition matrix \(k\) times, then take the weighted mean. The probabilities must sum to 1. For the matching uncertainty, see regimeSpread.
regime-modellingtime-series
def regimeMean(at, means) ≔ einsum("s,s->", at, means)
regimeSpread
regimeSpread(at, spreads) · at regime probabilities vector s · spreads standard deviation in each regime vector s → weighted variance scalar
Let \(p \in \Delta^{s-1}\) (the argument at) be the regime probabilities and \(\sigma \in \mathbb{R}^s\) (the argument spreads) the regime standard deviations. Then
The probability-weighted average of the regime variances. It does not add the dispersion of the regime means, so it equals the mixture variance only when the regimes share a mean.
Intuition The regime variances averaged with the regime probabilities as weights: the uncertainty within regimes. The full mixture variance also counts how far apart the regime means are.
Use it when you need an expected risk level from regime probabilities and the regimes have similar means, as volatility regimes of returns usually do. If the means differ, add \(\sum_k p_k(\mu_k - \bar\mu)^2\) yourself, with \(\bar\mu\) from regimeMean, or the risk will be understated.
regime-modellingrisk-measurement
def regimeSpread(at, spreads) ≔ einsum("s,s->", at, spreads · spreads)
regimeDrawn
regimeDrawn(at, means, spreads, z) · at regime probabilities vector s · means regime means vector s · spreads regime standard deviations vector s · z a standard normal draw scalar → an observation scalar
Let \(p \in \Delta^{s-1}\) (the argument at), \(\mu \in \mathbb{R}^s\) (the argument means) and \(\sigma \in \mathbb{R}^s\) (the argument spreads), and let \(z \in \mathbb{R}\) be a standard normal draw. Then
A draw with the weighted mean and weighted variance above. This is a single-normal approximation; to sample the mixture exactly, choose the regime first and use its own mean and spread.
Intuition The mixture is collapsed into a single normal with the blended mean and the within-regime variance, and then one draw is scaled. It is cheap and smooth, but it loses the mixture’s shape: two-humped or fat-tailed outcomes become a single bell curve.
Use it when you need quick approximate scenarios and the regimes are alike, or a smooth draw inside a larger model. For tail risk or realistic switching paths, sample exactly instead: pick the regime from at with a uniform draw (landedOn in Chains), use that regime’s own mean and spread, and move the regime along with the transition matrix.
regime-modellingsimulation
def regimeDrawn(at, means, spreads, z) ≔ regimeMean(at, means) + sqrt(max(regimeSpread(at, spreads), 0)) · z
Autoregression
The first-order autoregression \(x_{t+1} = \phi x_t + \sigma z_t\) is stationary when \(\lvert \phi \rvert < 1\) (for one lag this is the same as the characteristic root lying outside the unit circle). Its autocorrelation at lag \(k\) is \(\phi^k\) — geometric decay, the signature by which the process is recognised.
arStationary
arStationary(phi) · phi coefficient scalar → 1 or 0 scalar
Let \(\phi \in \mathbb{R}\) be the coefficient. Then
\[\operatorname{arStationary}(\phi) = [\,\lvert \phi \rvert < 1\,] .\]Intuition Each step keeps a share \(\phi\) of the current value. With \(\lvert\phi\rvert < 1\) old shocks shrink geometrically and the series has a stable level. At \(\lvert\phi\rvert = 1\) it is a random walk and every shock is permanent, and above that it explodes.
Use it when you check a fitted AR(1) before using arVariance, arHalfLife or long-horizon forecasts, all of which assume stationarity. An estimate just below 1 on a short series is hard to tell apart from a random walk, so back it with a unit-root test rather than this flag alone.
stability-checktime-series
def arStationary(phi) ≔ abs(phi) < 1
arVariance
arVariance(phi, sigma) · phi coefficient scalar · sigma noise standard deviation scalar → stationary variance scalar
Let \(\phi \in (-1, 1)\) be the coefficient and \(\sigma \ge 0\) the noise standard deviation. Then
\[\operatorname{arVariance}(\phi, \sigma) = \frac{\sigma^2}{\operatorname{nonzero}(1 - \phi^2)} .\]The stationary variance.
Intuition New noise adds \(\sigma^2\) each step while carrying over shrinks the existing variance by \(\phi^2\). At balance \(v = \phi^2 v + \sigma^2\), so the closer \(\lvert\phi\rvert\) is to 1, the larger the long-run swings from the same noise.
Use it when you turn innovation noise into the typical range of the level (a band of two standard deviations), start a simulation in equilibrium by drawing \(x_0\) with this variance, or check a simulated series. It needs \(\lvert\phi\rvert < 1\). The process here has mean 0, so add your own level. The continuous-time analogue is ouStationary in Diffusions.
time-seriesuncertainty-quantification
def arVariance(phi, sigma) ≔ sigma · sigma / nonzero(1 − phi · phi)
arCorrelation
arCorrelation(phi, k) · phi coefficient scalar · k lag scalar or vector → autocorrelation same shape as k
Let \(\phi \in (-1, 1)\) be the coefficient and \(k \in \mathbb{N}\) the lag (or a vector of lags, taken entrywise). Then
\[\operatorname{arCorrelation}(\phi, k) = \phi^{k} .\]The autocorrelation at lag \(k\).
Intuition A shock still has a share \(\phi^k\) of its size after \(k\) steps, so two values \(k\) apart share that much. A negative \(\phi\) gives correlations that alternate in sign.
Use it when you recognise an AR(1) in a sample autocorrelation plot (geometric decay), forecast several steps ahead (the \(k\)-step forecast of a zero-mean series is \(\phi^k x_t\)), or choose a spacing that gives nearly independent samples. Pass a vector of lags to get the whole curve. For an Ornstein–Uhlenbeck process sampled every \(\Delta t\), \(\phi\) is ouDecay(theta, dt) in Diffusions.
time-seriesparameter-estimation
def arCorrelation(phi, k) ≔ phi ^ k
arNext
arNext(x, phi, sigma, z) · x current value scalar · phi coefficient scalar · sigma noise standard deviation scalar · z a standard normal draw scalar → next value scalar
Let \(x_t \in \mathbb{R}\) (the argument x) be the current value, \(\phi \in \mathbb{R}\) the coefficient, \(\sigma \ge 0\) the noise standard deviation and \(z \in \mathbb{R}\) a standard normal draw. Then
One step of the process.
Intuition Keep a share \(\phi\) of today’s value and add a fresh shock: the simplest model in which the past matters. With \(0 < \phi < 1\) the series drifts back toward 0, and with a negative \(\phi\) it zigzags.
Use it when you simulate autocorrelated noise: sensor drift, demand deviations, error terms for testing a regression, or a baseline to compare volatility models against. Draw a fresh z each step. The series has mean 0, so add a constant for a level. If the noise size should itself cluster over time, take \(\sigma\) each step from garchNext instead of keeping it fixed.
simulationtime-series
def arNext(x, phi, sigma, z) ≔ phi · x + sigma · z
arHalfLife
arHalfLife(phi) · phi coefficient scalar → half-life in periods scalar
Let \(\phi \in (-1, 1)\) with \(\phi \ne 0\) be the coefficient. Then
\[\operatorname{arHalfLife}(\phi) = \frac{\ln 2}{\operatorname{nonzero}\big(-\ln \operatorname{nonzero}(\lvert \phi \rvert)\big)} .\]The number of periods for a deviation to halve.
Intuition The number of periods until an expected deviation shrinks to half, solving \(\lvert\phi\rvert^n = \tfrac12\). It uses \(\lvert\phi\rvert\), so for a negative \(\phi\) it measures how fast the zigzag dies down.
Use it when you interpret a fitted AR(1): how many days a price deviation or a spread takes to close halfway, or how long a demand shock echoes. It needs \(0 < \lvert\phi\rvert < 1\). At \(\phi = 0\) the guard gives about 0.001 periods, meaning no memory, and at \(\lvert\phi\rvert = 1\) an astronomically large value. The continuous-time version is halfLife in Diffusions.
mean-reversionmodel-interpretationtime-series
def arHalfLife(phi) ≔ ln(2) / nonzero(0 − ln(nonzero(abs(phi))))
Examples
The examples below live in one file, Clustering.hb, so the class is Clustering.
Filter a GARCH variance and read its long-run behaviour
With \(\omega = 0.1\), \(\alpha = 0.1\), \(\beta = 0.8\): a shock of 2 on a variance of 1 gives \(0.1 + 0.4 + 0.8 = 1.3\); the long-run variance is \(0.1/0.1 = 1\); and the kurtosis is \(3 \cdot 0.19/0.17\). A second model with persistence 0.95 has a shock half-life of about 13.5 periods.
open Prelude.Volatility
fn nextVariance(variance, shock, omega, alpha, beta) ↦ garchNext(variance, shock, omega, alpha, beta)
fn longRunVariance(omega, alpha, beta) ↦ unconditional(omega, alpha, beta)
fn halfLifeOf(alpha, beta) ↦ shockHalfLife(alpha, beta)
fn kurtosisOf(alpha, beta) ↦ garchKurtosis(alpha, beta)
fn isStationary(alpha, beta) ↦ garchStationary(alpha, beta)
using Examples;
var next = Clustering.NextVariance(1.0, 2.0, 0.1, 0.1, 0.8);
var longRun = Clustering.LongRunVariance(0.1, 0.1, 0.8);
var halfLife = Clustering.HalfLifeOf(0.05, 0.9);
var kurtosis = Clustering.KurtosisOf(0.1, 0.8);
var stable = Clustering.IsStationary(0.1, 0.8);
var integrated = Clustering.IsStationary(0.2, 0.8);
Output (verified, rounded to six digits):
garchNext = 1.3
unconditional = 1
halfLife = 13.513407
kurtosis = 3.352941
stationary = 1 / 0
The leverage effect in EGARCH
Starting from log variance 0 (variance 1), a shock of +1 and a shock of −1 with \(\gamma = -0.05\): the fall raises the next log variance by 0.1 more than the rise.
open Prelude.Volatility
fn nextLogVariance(logVariance, shock, omega, alpha, gamma, beta) ↦ egarchNext(logVariance, shock, omega, alpha, gamma, beta)
using Examples;
var afterRise = Clustering.NextLogVariance(0.0, 1.0, 0.0, 0.1, -0.05, 0.9);
var afterFall = Clustering.NextLogVariance(0.0, -1.0, 0.0, 0.1, -0.05, 0.9);
Output (verified, rounded to six digits):
egarch = -0.029788 / 0.070212
Correlated drivers and one Heston variance step
Two independent draws of 1 with \(\rho = 0.6\) give \(0.6 + 0.8 = 1.4\). A variance already at its long-run level of 0.04 moves only through the noise term: \(0.04 + 0.3 \cdot 0.2 \cdot 0.1 = 0.046\).
open Prelude.Volatility
fn correlatedDraw(first, second, rho) ↦ tied(first, second, rho)
fn hestonStep(v, kappa, theta, xi, dt, dw) ↦ hestonVariance(v, kappa, theta, xi, dt, dw)
using Examples;
var z2 = Clustering.CorrelatedDraw(1.0, 1.0, 0.6);
var v1 = Clustering.HestonStep(0.04, 2.0, 0.04, 0.3, 0.01, 0.1);
Output (verified):
tied = 1.4
heston = 0.046
How long a market regime lasts
A calm regime that stays put with probability 0.9 lasts ten periods on average; a turbulent one at 0.8 lasts five. With regime probabilities 0.75 and 0.25 and means 0.02 and −0.04, the mixture mean is 0.005. The transition matrix goes in as a Tensor because the result is a vector.
open Prelude.Volatility
fn regimeLengths(transitions) ↦ regimeDuration(transitions)
fn mixedMean(at, means) ↦ regimeMean(at, means)
using MLambda.Hilbert.Runtime;
using Examples;
var durations = Clustering.RegimeLengths((Tensor)new double[,] { { 0.9, 0.1 }, { 0.2, 0.8 } });
var mean = Clustering.MixedMean(new[] { 0.75, 0.25 }, new[] { 0.02, -0.04 });
Output (verified):
durations = [10, 5]
regimeMean = 0.005
Variance and half-life of an AR(1)
At \(\phi = 0.5\) and \(\sigma = 1\) the stationary variance is \(1/0.75\) and a deviation halves every period.
open Prelude.Volatility
fn arSpread(phi, sigma) ↦ arVariance(phi, sigma)
fn arLife(phi) ↦ arHalfLife(phi)
using Examples;
var variance = Clustering.ArSpread(0.5, 1.0);
var halfLife = Clustering.ArLife(0.5);
Output (verified, rounded to six digits):
arVariance = 1.333333
arHalfLife = 1