Diffusions
Stochastic differential equations — Euler–Maruyama and Milstein steps, with the closed forms of Brownian motion, GBM, Ornstein–Uhlenbeck, CIR, Merton jumps and fractional Brownian motion.
This module covers processes in continuous time with a continuous state, written as a stochastic differential equation \(dX = a(X)\,dt + b(X)\,dW\). Two discretisation schemes serve the whole family: Euler–Maruyama and Milstein. Brownian motion, geometric Brownian motion, Ornstein–Uhlenbeck, the square-root (CIR) diffusion and Merton’s jump diffusion differ only in the drift \(a\) and diffusion \(b\) they pass to the scheme, and each comes with its drift and diffusion terms and its exact moments. Where a closed-form solution exists it sits beside the scheme, so a simulated path can be checked against the process it approximates.
Conventions. A step function takes the current state x, the drift a and diffusion b evaluated at x, the time step dt, and the Brownian increment dw. Build dw from a standard normal draw with increment(z, dt), which is \(z\sqrt{dt}\); the normal draws come from the caller or, in a process declaration, from a gaussian stream (see the Brownian, Geometric, Ornstein, Square, Merton, Fractional, Maruyama and Milstein sample models). Euler–Maruyama converges with strong order ½ and weak order 1; Milstein has strong order 1.
Limitations. Fractional Brownian motion is generated exactly through a Cholesky factor of the n×n covariance, so memory grows as \(n^2\) and time faster; it suits paths of hundreds to a few thousand steps. The Ornstein–Uhlenbeck fit is a least-squares line whose averages are taken over the whole series rather than only over the paired points, which biases it on short series; with a few hundred observations it agrees closely with an exact regression (see the example). quadraticVariation builds an n×n mask.
In the statements below, \(\operatorname{nonzero}(v) = v + [\,\lvert v \rvert < 10^{-300}\,] \cdot 10^{-300}\) (applied entrywise) is the divisor guard from Chains; it moves only entries of magnitude below \(10^{-300}\), in practice exact zeros. Where a definition reads positions in a series, positions are counted from zero.
open Prelude.Diffusions
The two schemes
Euler–Maruyama and Milstein advance the state over one step of length \(\Delta t\):
\[X_{k+1} = X_k + a\,\Delta t + b\,\Delta W, \qquad X_{k+1}^{\text{Milstein}} = X_{k+1}^{\text{Euler}} + \tfrac12\, b\, b' \big(\Delta W^2 - \Delta t\big)\]where \(b'\) is the derivative of the diffusion coefficient with respect to the state.
euler
euler(x, a, b, dt, dw) · x state · a drift at x · b diffusion at x · dt time step · dw Brownian increment → next state
Let \(x \in \mathbb{R}\) be the state, \(a = a(x)\) and \(b = b(x)\) the drift and diffusion evaluated there, \(\Delta t > 0\) (the argument dt) the time step and \(\Delta W \in \mathbb{R}\) (the argument dw) the Brownian increment. Then
One Euler–Maruyama step. Works elementwise on tensors of states.
Intuition Over a short step the state is pushed by the drift for the elapsed time and jostled by a random kick whose size grows like \(\sqrt{\Delta t}\). The drift and diffusion are frozen at their values at the start of the step, which is accurate when they barely change during it.
Use it when you simulate any SDE path: pass the model’s drift and diffusion at the current state (gbmDrift, ouDrift, cirDiffusion, …) and a draw from increment. It is the default for Monte Carlo prices and risk, where only averages matter (weak order 1). Single paths converge only like \(\sqrt{\Delta t}\), so when the path matters and the noise depends on the state use milstein. Where a closed form exists (gbmExact) it beats any step size.
simulationasset-pricinguncertainty-quantification
def euler(x, a, b, dt, dw) ≔ x + a · dt + b · dw
milstein
milstein(x, a, b, db, dt, dw) · x state · a drift at x · b diffusion at x · db derivative of the diffusion at x · dt time step · dw Brownian increment → next state
Let \(x \in \mathbb{R}\) be the state, \(a = a(x)\), \(b = b(x)\) and \(b' = \partial b / \partial x\) (the argument db) evaluated there, \(\Delta t > 0\) the time step and \(\Delta W \in \mathbb{R}\) the Brownian increment. Then
One Milstein step: the Euler step plus a correction. When the diffusion does not depend on the state (\(b' = 0\)) it equals the Euler step.
Intuition Euler treats the noise level as constant over the step. When it depends on the state, the kick itself changes the noise level during the step. Milstein adds the next Itô–Taylor term for that. The term has mean zero, so averages are unchanged, but single paths track the true path more closely.
Use it when path accuracy matters and the diffusion depends on the state (GBM with \(b' = \sigma\), local volatility): path-dependent payoffs, barrier checks, coupling fine and coarse paths in multilevel Monte Carlo. You supply \(b'\) yourself. With additive noise (Ornstein–Uhlenbeck, Brownian motion) \(b' = 0\) and it is just euler. Avoid it for the square-root diffusion near zero, where \(b'\) blows up.
simulationoption-pricing
def milstein(x, a, b, db, dt, dw) ≔ euler(x, a, b, dt, dw) + 0.5 · b · db · (dw · dw − dt)
increment
increment(z, dt) · z standard normal draw · dt time step → Brownian increment
Let \(z \in \mathbb{R}\) be a standard normal draw and \(\Delta t \ge 0\) (the argument dt) the time step. Then
The Brownian increment, with variance \(\Delta t\).
Intuition Brownian motion’s change over \(\Delta t\) is normal with variance \(\Delta t\), so its size scales with \(\sqrt{\Delta t}\), not \(\Delta t\). That square root is why noise outweighs drift on short steps.
Use it when you turn standard normal draws into the dw argument of euler, milstein, mertonStep, or the Heston steps in Volatility. Use the same dt as the step, and a fresh independent draw for every step and path. For correlated drivers, correlate the normals first (tied in Volatility), then scale.
simulationuncertainty-quantification
def increment(z, dt) ≔ z · sqrt(dt)
Brownian motion
Brownian motion with volatility \(\sigma\) has variance \(\sigma^2 t\) that grows without limit. Its quadratic variation \(\sum (\Delta W)^2\) converges to \(t\) rather than to zero, which is why Itô calculus needs a second-order term. By the reflection principle the running maximum exceeds a level \(\ell > 0\) twice as often as the endpoint does, and with a drift \(\mu\) carrying the process away from the level it is ever reached with probability \(e^{-2\mu\ell/\sigma^2}\):
\[P\Big(\max_{s \le t} W_s \ge \ell\Big) = 2\,P(W_t \ge \ell), \qquad P(\text{ever reach } \ell) = \min\big(e^{-2\mu\ell/\sigma^2},\, 1\big)\]brownianSpread
brownianSpread(sigma, t) · sigma volatility σ · t time → variance scalar
Let \(\sigma \ge 0\) (the argument sigma) be the volatility and \(t \ge 0\) a time. Then
The variance of Brownian motion at time \(t\).
Intuition Independent kicks add their variances, so uncertainty builds up linearly in time and the standard deviation grows like \(\sigma\sqrt{t}\). Nothing pulls the process back, so the spread never stops growing.
Use it when you need to know how uncertain a random walk is at a horizon: the error of a sensor that drifts randomly, the spread of a log-price over \(t\), or the square-root-of-time rule for scaling daily volatility to annual. It is also a check on a simulated ensemble. If the quantity is pulled back toward a level its variance levels off, so use ouSpread.
uncertainty-quantificationrisk-measurementmodel-validation
def brownianSpread(sigma, t) ≔ sigma · sigma · t
quadraticVariation
quadraticVariation(path) · path sampled values vector n → quadratic variation scalar
Let \(x = (x_0, \dots, x_{n-1}) \in \mathbb{R}^n\) (the argument path) be sampled values, positions counted from zero, and let \(x_{-1} = 0\). Then
where \(x_{i-1}\) is read from \(\operatorname{previousOf}(x)\) in Paths.
The sum of squared increments. The value before the first sample is taken as 0, so a path that starts at 0 (as a process weight does) is measured from its origin; for a path that starts elsewhere, the first term includes the square of its starting value.
Intuition Add up the squared moves along a path. For a smooth curve this goes to zero as sampling gets finer. For Brownian motion it converges to \(\sigma^2 t\), because each squared move has size about \(\Delta t\) and there are \(t/\Delta t\) of them. It measures how much the path jiggles, whatever its endpoint.
Use it when you estimate realised variance from high-frequency data (pass log-prices), check that a simulated path has the diffusion coefficient you meant, or tell rough noise from a smooth trend. Start the path at 0, or its starting value squared is added. Microstructure noise on very fine grids inflates it. It builds an n×n mask, so memory grows as \(n^2\).
parameter-estimationmodel-validation
def quadraticVariation(path) ≔ einsum("i->", (path − previousOf(path)) · (path − previousOf(path)))
reachesLevel
reachesLevel(level, sigma, t) · level ℓ > 0 · sigma volatility σ · t time → probability scalar
Let \(\ell > 0\) (the argument level), \(\sigma > 0\) and \(t > 0\). Then
with \(\Phi\) computed by normalCdf from Distributions.
The probability that driftless Brownian motion started at 0 reaches level by time \(t\).
Intuition The reflection principle: every path that touches \(\ell\) and ends below it has a mirror twin, flipped at the first touch, that ends above it and is just as likely. So “touched \(\ell\) by \(t\)” is exactly twice as likely as “above \(\ell\) at \(t\)”, which is a normal tail.
Use it when you need the chance of crossing a threshold at any time during a period, not just at the end: a stop-loss or knock-out barrier, a stock-out, a limit being exceeded. It assumes no drift, a start at 0 (use the distance to the level) and continuous monitoring. A simulation checked only at grid points reports fewer hits. With a drift away from the level and no time limit, use everReaches.
barrier-crossingrisk-measurementoption-pricing
def reachesLevel(level, sigma, t) ≔ 2 · (1 − normalCdf(0, sigma · sqrt(t), level))
everReaches
everReaches(level, drift, sigma) · level ℓ > 0 · drift speed away from the level μ · sigma volatility σ → probability scalar
Let \(\ell > 0\) (the argument level), \(\mu \in \mathbb{R}\) (the argument drift, positive when it carries the process away from the level) and \(\sigma > 0\). Then
The probability that Brownian motion ever reaches level when its drift carries it away from the level at speed drift. A drift of zero or toward the level gives 1.
Intuition A drift away from the level wins in the long run, but early noise can carry the process across before the drift takes over. The chance shrinks exponentially as the level gets farther away, the drift gets stronger or the noise gets quieter.
Use it when you estimate ruin or exceedance with no time limit: an insurer’s surplus drifting up by its premium margin ever falling to a loss level, a reservoir ever overflowing, a workload ever reaching a bound. Measure the level as a distance from the start. A drift of zero or toward the level returns 1, meaning certain. For a finite horizon with no drift, use reachesLevel.
barrier-crossingrisk-measurementinsurance-claims
def everReaches(level, drift, sigma) ≔ min(exp(0 − 2 · drift · level / nonzero(sigma · sigma)), 1)
Geometric Brownian motion
\(dX = \mu X\,dt + \sigma X\,dW\) has the exact solution
\[X_t = X_0 \exp\Big(\big(\mu - \tfrac{\sigma^2}{2}\big)t + \sigma\sqrt{t}\,Z\Big)\]The mean grows at rate \(\mu\) but the median at \(\mu - \sigma^2/2\), so a process with a positive expected return can still have a falling median: that is what compounding does to variance.
gbmDrift
gbmDrift(x, mu) · x state · mu growth rate μ → drift
Let \(x \in \mathbb{R}\) be the state and \(\mu \in \mathbb{R}\) (the argument mu) the growth rate. Then
The GBM drift, to pass as a to a scheme.
Intuition Growth in proportion to size: a price, balance or population earning rate \(\mu\) grows by \(\mu x\) per unit time. What stays constant is the percentage change, not the absolute change.
Use it when you pass a to euler or milstein to simulate asset prices, compounding balances or unchecked population growth, where the quantity is positive and its moves scale with its level. For large steps gbmExact has no discretisation error. For a quantity pulled toward a level, use ouDrift.
asset-pricingsimulation
def gbmDrift(x, mu) ≔ mu · x
gbmDiffusion
gbmDiffusion(x, sigma) · x state · sigma volatility σ → diffusion
Let \(x \in \mathbb{R}\) be the state and \(\sigma \ge 0\) (the argument sigma) the volatility. Then
The GBM diffusion, to pass as b to a scheme (its derivative, for Milstein, is \(\sigma\)).
Intuition Noise in proportion to size: a stock at 200 moves twice as much as one at 100 on the same percentage shock. That makes \(\ln x\) a plain Brownian motion with drift, and in continuous time \(x\) never reaches zero.
Use it when you pass b to a scheme for price or growth paths with constant percentage volatility, as in Black–Scholes. For milstein pass db = sigma. A discrete Euler step can still go negative when \(\sigma\sqrt{\Delta t}\) is large, and gbmExact avoids that. If the volatility itself should move, use the Heston steps in Volatility.
asset-pricingoption-pricingsimulation
def gbmDiffusion(x, sigma) ≔ sigma · x
gbmExact
gbmExact(x0, mu, sigma, t, z) · x0 starting value · mu μ · sigma σ · t time · z standard normal draw → value at time t
Let \(X_0 > 0\) (the argument x0), \(\mu \in \mathbb{R}\), \(\sigma \ge 0\), \(t \ge 0\) and let \(z \in \mathbb{R}\) be a standard normal draw. Then
A draw of \(X_t\) from the exact solution, using the single normal z.
Intuition The log of the price is Brownian motion with drift \(\mu - \sigma^2/2\), so the price at \(t\) is the start times the exponential of one normal variable. A single normal draw gives an exact sample at any horizon, with no stepping error and never a negative value.
Use it when you price European options by Monte Carlo, generate terminal price scenarios for value-at-risk, or check where Euler paths end up. For a whole path, apply it grid step by grid step with fresh draws: that is still exact at every grid point. It needs \(X_0 > 0\). Models with no closed form (stochastic volatility, mertonStep jumps) must be stepped with a scheme.
option-pricingasset-pricingrisk-measurement
def gbmExact(x0, mu, sigma, t, z) ≔ x0 · exp((mu − sigma · sigma / 2) · t + sigma · sqrt(t) · z)
gbmMean
gbmMean(x0, mu, t) · x0 starting value · mu μ · t time → mean scalar
Let \(X_0 \in \mathbb{R}\) (the argument x0), \(\mu \in \mathbb{R}\) and \(t \ge 0\). Then
The expected value.
Intuition Averaged over all paths, value compounds at exactly \(\mu\). The lognormal’s long upper tail makes up for the \(-\sigma^2/2\) in the exponent, so the mean is held up by a few very lucky paths.
Use it when you need an expected value: a forward price, an expected portfolio value, or the target a Monte Carlo average should hit. Do not read it as the typical outcome. With high volatility or a long horizon most paths end well below it, so use gbmMedian for the typical path and gbmSpread for the dispersion.
asset-pricingmodel-validation
def gbmMean(x0, mu, t) ≔ x0 · exp(mu · t)
gbmMedian
gbmMedian(x0, mu, sigma, t) · x0 starting value · mu μ · sigma σ · t time → median scalar
Let \(X_0 > 0\) (the argument x0), \(\mu \in \mathbb{R}\), \(\sigma \ge 0\) and \(t \ge 0\). Then
The median of \(X_t\).
Intuition Half of all paths end below this value. It grows at \(\mu - \sigma^2/2\) because volatility drags on compounding: up 10% and then down 10% leaves you at 99%. When \(\sigma^2 > 2\mu\) the median falls even while the mean rises.
Use it when you want the outcome a single holder or a single realisation should expect: long-horizon wealth projections, a “typical” scenario, or explaining why a volatile asset with a positive expected return loses money for most holders. The gap to gbmMean shows how much the mean depends on rare large outcomes.
asset-pricingrisk-measurement
def gbmMedian(x0, mu, sigma, t) ≔ x0 · exp((mu − sigma · sigma / 2) · t)
gbmSpread
gbmSpread(x0, mu, sigma, t) · x0 starting value · mu μ · sigma σ · t time → variance scalar
Let \(X_0 \in \mathbb{R}\) (the argument x0), \(\mu \in \mathbb{R}\), \(\sigma \ge 0\) and \(t \ge 0\). Then
The variance of \(X_t\).
Intuition The variance of a lognormal. It grows like \(e^{\sigma^2 t} - 1\) relative to the squared mean, because compounding has no ceiling on the upside while the downside stops at zero.
Use it when you need the dispersion of a future value for a risk budget, a band on a price forecast, or a check on simulation output. The distribution is strongly skewed, so the mean plus or minus two standard deviations is a poor interval and its lower end can go negative. For quantiles, use the normal law of \(\ln X_t\): gbmExact with z set to the normal quantile you want.
risk-measurementuncertainty-quantificationmodel-validation
def gbmSpread(x0, mu, sigma, t) ≔ x0 · x0 · exp(2 · mu · t) · (exp(sigma · sigma · t) − 1)
Ornstein–Uhlenbeck
\(dX = \theta(\mu - X)\,dt + \sigma\,dW\) pulls the state back toward \(\mu\) in proportion to how far it has strayed. Unlike Brownian motion it has a stationary distribution, \(N(\mu, \sigma^2/2\theta)\), and its moments from a fixed start are
\[E[X_t] = \mu + (x_0 - \mu)e^{-\theta t}, \qquad \operatorname{Var} X_t = \frac{\sigma^2}{2\theta}\big(1 - e^{-2\theta t}\big)\]ouDrift
ouDrift(x, theta, mu) · x state · theta reversion speed θ · mu long-run mean μ → drift
Let \(x \in \mathbb{R}\) be the state, \(\theta > 0\) (the argument theta) the reversion speed and \(\mu \in \mathbb{R}\) (the argument mu) the long-run mean. Then
The mean-reverting drift, to pass as a to a scheme.
Intuition A spring toward \(\mu\): the farther the state strays, the harder it is pulled back, and \(\theta\) is the stiffness. With constant noise on top, the state wanders around \(\mu\) instead of drifting off.
Use it when you simulate a quantity with a natural level: interest rates (Vasicek), spreads in pairs trading, temperature around a seasonal mean, a velocity with friction. Pass it as a with a constant b = sigma. Keep \(\theta\,\Delta t\) well below 1, because Euler overshoots the mean above 1 and diverges above 2. If the quantity must stay positive, pair it with cirDiffusion.
mean-reversioninterest-rate-modellingsimulation
def ouDrift(x, theta, mu) ≔ theta · (mu − x)
ouMean
ouMean(theta, mu, x0, t) · theta θ · mu μ · x0 starting value · t time → mean scalar
Let \(\theta > 0\), \(\mu \in \mathbb{R}\), \(x_0 \in \mathbb{R}\) (the argument x0) and \(t \ge 0\). Then
The expected value at time \(t\).
Intuition The expected gap to \(\mu\) shrinks by the same fraction per unit time, an exponential approach like a hot cup cooling to room temperature. The noise has no effect on the average.
Use it when you forecast a mean-reverting variable, such as the expected rate, spread or temperature at a horizon, build the expected path of a scenario, or check simulated averages. For the uncertainty around the forecast, use ouSpread. The square-root diffusion has the same mean, and cirMean is this function under that name.
mean-reversioninterest-rate-modellingmodel-validation
def ouMean(theta, mu, x0, t) ≔ mu + (x0 − mu) · exp(0 − theta · t)
ouSpread
ouSpread(theta, sigma, t) · theta θ · sigma σ · t time → variance scalar
Let \(\theta > 0\), \(\sigma \ge 0\) and \(t \ge 0\). Then
\[\operatorname{ouSpread}(\theta, \sigma, t) = \frac{\sigma^2}{\operatorname{nonzero}(2\theta)} \big( 1 - e^{-2\theta t} \big) .\]The variance at time \(t\) from a fixed start.
Intuition Noise pumps variance in at rate \(\sigma^2\) while the pull toward \(\mu\) drains it at rate \(2\theta\). From a known start, the variance climbs and levels off where the two balance, getting most of the way there within a few multiples of \(1/(2\theta)\).
Use it when you need a confidence band on a mean-reverting forecast (the state at \(t\) is normal with ouMean and this variance), an exact transition for simulation or a likelihood (use ouMean(theta, mu, x, dt) and ouSpread(theta, sigma, dt) for one step of length dt), or to see how soon a forecast stops being informative. Long horizons converge to ouStationary.
mean-reversionuncertainty-quantificationlikelihood-evaluation
def ouSpread(theta, sigma, t) ≔ sigma · sigma / nonzero(2 · theta) · (1 − exp(0 − 2 · theta · t))
ouStationary
ouStationary(theta, sigma) · theta θ · sigma σ → variance scalar
Let \(\theta > 0\) and \(\sigma \ge 0\). Then
\[\operatorname{ouStationary}(\theta, \sigma) = \frac{\sigma^2}{\operatorname{nonzero}(2\theta)} .\]The stationary variance.
Intuition The long-run balance between noise pumping variance in (\(\sigma^2\)) and the restoring pull draining it (\(2\theta\)). A stiffer spring or quieter noise gives a tighter band around \(\mu\).
Use it when you want the band a mean-reverting variable normally lives in (\(\mu\) plus or minus two standard deviations), to start a simulation in equilibrium by drawing \(x_0\) from \(N(\mu, \sigma^2/2\theta)\) instead of fixing it, or to recover \(\sigma\) from an observed long-run standard deviation and a known \(\theta\). As \(\theta\) nears 0 it blows up, because the process is then nearly Brownian motion and has no long-run band.
mean-reversionrisk-measurementparameter-estimation
def ouStationary(theta, sigma) ≔ sigma · sigma / nonzero(2 · theta)
halfLife
halfLife(theta) · theta θ → time scalar
Let \(\theta > 0\) be the reversion speed. Then
\[\operatorname{halfLife}(\theta) = \frac{\ln 2}{\operatorname{nonzero}(\theta)} .\]The time for a deviation from the mean to halve on average.
Intuition \(\theta\) is an abstract rate. The half-life turns it into a duration: how long it takes, on average, for a deviation from \(\mu\) to shrink by half. As with radioactive decay, it is the same whatever the size of the starting gap.
Use it when you choose or explain \(\theta\) in practical terms: how many days a pairs-trading spread takes to close halfway, how quickly rates revert, whether reversion is fast enough to trade within a holding period. The answer is in the time unit of \(\theta\). The discrete-time counterparts are shockHalfLife (GARCH) and arHalfLife (AR(1)) in Volatility.
mean-reversionmodel-interpretation
def halfLife(theta) ≔ ln(2) / nonzero(theta)
ouDecay
ouDecay(theta, t) · theta θ · t lag → autocorrelation scalar
Let \(\theta > 0\) and \(t \ge 0\) a lag. Then
\[\operatorname{ouDecay}(\theta, t) = e^{-\theta t} .\]The stationary autocorrelation at lag \(t\).
Intuition In equilibrium, how much today’s deviation still predicts the deviation a lag \(t\) later. The memory fades at rate \(\theta\), so values more than a few multiples of \(1/\theta\) apart are nearly independent.
Use it when you match a model to data through the autocorrelation (the log of the sample autocorrelation against the lag has slope \(-\theta\)), choose a sampling interval that gives nearly independent observations, or convert to discrete time: sampled every \(\Delta t\), the process is an AR(1) with coefficient ouDecay(theta, dt) (see arCorrelation). To fit \(\theta\) by regression from a series, use thetaFrom.
mean-reversionparameter-estimationtime-series
def ouDecay(theta, t) ≔ exp(0 − theta · t)
The square-root diffusion
The Cox–Ingersoll–Ross process \(dX = \theta(\mu - X)\,dt + \sigma\sqrt{X}\,dW\) has the same mean as Ornstein–Uhlenbeck, but its noise vanishes at zero. Feller’s condition \(2\theta\mu \ge \sigma^2\) keeps it strictly positive. Below that, the process can touch zero, a naive discretisation takes the square root of a negative number, and every later value is NaN. Full truncation, which uses the positive part of the state inside the square root, is the standard repair.
cirDiffusion
cirDiffusion(x, sigma) · x state · sigma σ → diffusion
Let \(x \in \mathbb{R}\) be the state and \(\sigma \ge 0\). Then
\[\operatorname{cirDiffusion}(x, \sigma) = \sigma \sqrt{\max(x, 0)} .\]The square-root diffusion, with full truncation built in.
Intuition Noise proportional to the square root of the level: large when the level is high, fading to zero as it nears zero, so near zero the pull back toward \(\mu\) wins. The \(\max(x, 0)\) inside stops a discrete step that dips below zero from producing NaN.
Use it when you simulate a quantity that must stay non-negative and whose volatility grows with its level: short interest rates (the CIR model), a stochastic variance, a default intensity. Use ouDrift as a, evaluated at fullTruncation(x), to get the full-truncation scheme. Check fellerHolds first, because if it fails paths touch zero often and the step size matters more. Prefer Euler to Milstein here.
interest-rate-modellingmean-reversionnumerical-stability
def cirDiffusion(x, sigma) ≔ sigma · sqrt(max(x, 0))
fullTruncation
fullTruncation(x) · x state → positive part
Let \(x \in \mathbb{R}\) be the state. Then
\[\operatorname{fullTruncation}(x) = \max(x, 0) .\]The positive part, to use in place of the state wherever the scheme needs it non-negative (for example inside the drift).
Intuition The continuous process never goes below zero, but a discrete step can overshoot. Using the positive part wherever a formula needs a non-negative input lets the path recover, because the drift pushes it back up, instead of poisoning every later value with NaN.
Use it when you step a CIR rate or a Heston-type variance with Euler. Use it inside both evaluations, for example euler(x, ouDrift(fullTruncation(x), theta, mu), cirDiffusion(x, sigma), dt, dw). The stored state may stay slightly negative between steps, so apply it again before reporting a rate or variance. hestonVariance in Volatility instead clamps the result of each step at zero.
numerical-stabilityinterest-rate-modellingsimulation
def fullTruncation(x) ≔ max(x, 0)
fellerHolds
fellerHolds(theta, mu, sigma) · theta θ · mu μ · sigma σ → flag 1 or 0
Let \(\theta > 0\), \(\mu > 0\) and \(\sigma \ge 0\). Then
\[\operatorname{fellerHolds}(\theta, \mu, \sigma) = [\, 2\theta\mu \ge \sigma^2 \,] .\]1 when Feller’s condition holds, so the process stays strictly positive.
Intuition Near zero the drift pushes up at about \(\theta\mu\) while the noise shrinks with the level. If the push is strong enough relative to the noise, \(2\theta\mu \ge \sigma^2\), the process gets close to zero but never reaches it.
Use it when you validate calibrated CIR parameters before simulating or pricing. A 0 does not make the model invalid, since zero simply becomes reachable. But discretisation schemes lose accuracy near zero, so use smaller steps with fullTruncation and expect more probability mass near zero. The same test under Heston’s parameter names is hestonFeller in Volatility.
stability-checknumerical-stabilityinterest-rate-modelling
def fellerHolds(theta, mu, sigma) ≔ 2 · theta · mu ≥ sigma · sigma
cirMean
cirMean(theta, mu, x0, t) · theta θ · mu μ · x0 starting value · t time → mean scalar
Let \(\theta > 0\), \(\mu > 0\), \(x_0 \ge 0\) (the argument x0) and \(t \ge 0\). Then
The expected value at time \(t\), identical to the Ornstein–Uhlenbeck mean.
Intuition Square-root noise changes the shape of the distribution and keeps it positive, but not its average. The drift is still linear in \(x\), so the expected value approaches \(\mu\) exponentially, exactly as in Ornstein–Uhlenbeck.
Use it when you forecast the expected short rate or variance at a horizon under a square-root model, or check a truncated simulation’s average for bias. It is ouMean under a name that states the model. The variance is not the Ornstein–Uhlenbeck one: its long-run value is cirStationary.
interest-rate-modellingmean-reversionmodel-validation
def cirMean(theta, mu, x0, t) ≔ ouMean(theta, mu, x0, t)
cirStationary
cirStationary(theta, mu, sigma) · theta θ · mu μ · sigma σ → variance scalar
Let \(\theta > 0\), \(\mu > 0\) and \(\sigma \ge 0\). Then
\[\operatorname{cirStationary}(\theta, \mu, \sigma) = \frac{\sigma^2 \mu}{\operatorname{nonzero}(2\theta)} .\]The stationary variance.
Intuition The noise variance at the long-run level is \(\sigma^2\mu\) rather than \(\sigma^2\), so the balance with the pull \(2\theta\) gains a factor \(\mu\) compared with ouStationary. A higher long-run level means wider fluctuations.
Use it when you need the typical band of a CIR rate or a stochastic variance around its mean, want to calibrate \(\sigma\) from an observed long-run variance, or are checking a long simulation. The stationary law is a gamma distribution skewed to the right, so a symmetric two-standard-deviation band is only rough. It needs \(\theta > 0\).
interest-rate-modellingmean-reversionparameter-estimation
def cirStationary(theta, mu, sigma) ≔ sigma · sigma · mu / nonzero(2 · theta)
Jumps on top of a diffusion
Merton’s model adds a compound Poisson process of proportional jumps to geometric Brownian motion. A pure diffusion has kurtosis 3 whatever its parameters, while real returns have fatter tails; jumps of standard deviation \(s\) arriving at rate \(\lambda\) raise the kurtosis of the log-return over time \(t\) to
\[3 + \frac{3\lambda t\, s^4}{\big(\sigma^2 t + \lambda t\, s^2\big)^2}\]jumpArrived
jumpArrived(u, rate, dt) · u uniform draw · rate jump rate λ · dt time step → flag 1 or 0
Let \(u \in [0, 1)\) be a uniform draw, \(\lambda \ge 0\) (the argument rate) the jump rate and \(\Delta t > 0\) (the argument dt) the time step. Then
1 when a jump occurs in this step, with probability \(\lambda\,\Delta t\). This is the small-step approximation, so keep \(\lambda\,\Delta t\) well below one.
Intuition In a short step a Poisson process fires with probability of about \(\lambda\,\Delta t\) and almost never twice. Comparing a uniform draw with \(\lambda\,\Delta t\) is a coin flip with exactly that bias.
Use it when you add sudden events to a time-stepped simulation: price gaps (with mertonStep), outages, defaults, catastrophe losses. It allows at most one event per step and undercounts once \(\lambda\,\Delta t\) is not small (the exact chance of at least one event is \(1 - e^{-\lambda\Delta t}\)), so keep \(\lambda\,\Delta t\) to about 0.1 or less. For exact event times, use waitFor in Jumps.
sudden-shockssimulation
def jumpArrived(u, rate, dt) ≔ u < rate · dt
mertonStep
mertonStep(x, mu, sigma, dt, dw, jumped, size) · x state · mu μ · sigma σ · dt time step · dw Brownian increment · jumped jump flag 1 or 0 · size proportional jump size → next state
Let \(x \in \mathbb{R}\) be the state, \(\mu \in \mathbb{R}\), \(\sigma \ge 0\), \(\Delta t > 0\) (the argument dt), \(\Delta W \in \mathbb{R}\) (the argument dw) the Brownian increment, \(\eta \in \{0, 1\}\) (the argument jumped) the jump flag and \(y \in \mathbb{R}\) (the argument size) the proportional jump size. Then
One step of Merton’s jump diffusion: an Euler step of GBM plus \(x y\) when a jump arrived.
Intuition Ordinary GBM movement plus, now and then, a sudden proportional gap: when a jump arrives the price is multiplied by \(1 + y\). Most of the time the price diffuses, and occasionally it gaps, as on earnings news or in a crash.
Use it when you simulate assets with gap risk (equity crashes, devaluations, credit events) to price out-of-the-money options or run stress tests. Get jumped from jumpArrived and draw size from your jump distribution. The jumps are not compensated, so the expected growth rate is \(\mu + \lambda E[y]\). Pass \(\mu - \lambda E[y]\) if \(\mu\) should be the total expected return. A size below \(-1\) makes the price negative.
sudden-shocksoption-pricingasset-pricing
def mertonStep(x, mu, sigma, dt, dw, jumped, size) ≔ euler(x, gbmDrift(x, mu), gbmDiffusion(x, sigma), dt, dw) + x · jumped · size
mertonKurtosis
mertonKurtosis(sigma, rate, jumpSpread, t) · sigma σ · rate jump rate λ · jumpSpread jump standard deviation s · t time → kurtosis scalar
Let \(\sigma \ge 0\), \(\lambda \ge 0\) (the argument rate), \(s \ge 0\) (the argument jumpSpread) the standard deviation of a zero-mean normal jump, and \(t > 0\). Then
The kurtosis of the return over time \(t\) with zero-mean normal jumps; 3 when there are no jumps.
Intuition Kurtosis measures how much of the variance comes from rare large moves. A diffusion alone gives normal log-returns and kurtosis 3. Jumps add a few big moves, and the excess is largest when jumps make up a big share of the variance. The excess shrinks like \(1/t\), because over longer horizons many jumps average out toward a normal shape.
Use it when you calibrate jump rate and size to the fat tails seen in returns, or explain why daily returns look more fat-tailed than yearly ones. It assumes zero-mean normal jumps. If the fat tails might come from changing variance rather than from jumps, compare garchKurtosis in Volatility.
sudden-shocksrisk-measurementparameter-estimation
def mertonKurtosis(sigma, rate, jumpSpread, t) ≔ 3 + rate · t · 3 · jumpSpread ^ 4 / nonzero((sigma · sigma · t + rate · t · jumpSpread · jumpSpread) ^ 2)
Fractional Brownian motion
Fractional Brownian motion with Hurst exponent \(H \in (0, 1)\) has covariance
\[\operatorname{Cov}(B_s, B_t) = \tfrac12\big(s^{2H} + t^{2H} - |s - t|^{2H}\big)\]which depends on both times and not on their difference alone, so the process is not Markov: it remembers its history. \(H > \tfrac12\) is persistent (a rise tends to be followed by a rise), \(H < \tfrac12\) antipersistent, and \(H = \tfrac12\) is ordinary Brownian motion. A path is generated exactly as \(L z\), where \(L\) is the Cholesky factor of the covariance at times \(1, \dots, n\).
fbmCovariance
fbmCovariance(z, h) · z any vector whose length n sets the number of times vector n · h Hurst exponent H → covariance n×n
Let \(z \in \mathbb{R}^n\) be any vector (only its length \(n\) is used), \(H \in (0, 1)\) (the argument h) the Hurst exponent, and \(\tau_i = i + 1\) for positions \(i = 0, \dots, n-1\). Then
The covariance matrix of fractional Brownian motion at the integer times \(1, 2, \dots, n\).
Intuition Ordinary Brownian motion’s covariance is \(\min(s, t)\): two times share only their common history. With \(H \ne \tfrac12\) even increments far apart are correlated, positively for \(H > \tfrac12\) (rises follow rises) and negatively below (moves reverse). The value at time \(t\) has variance \(t^{2H}\).
Use it when you need the covariance matrix of a long-memory model: to build a path with fbmPath, to write a Gaussian likelihood for estimating \(H\), or to inspect the dependence. The times are \(1, \dots, n\), so for spacing \(\Delta t\) multiply by \(\Delta t^{2H}\). Only the length of z matters. At \(H = \tfrac12\) it reduces to \(\min(s, t)\), which makes a good sanity check.
long-memorylikelihood-evaluationtime-series
def fbmCovariance(z, h) ≔ 0.5 · (held(fbmTimes(z)) ^ (2 · h) + laid(fbmTimes(z)) ^ (2 · h) − abs(held(fbmTimes(z)) − laid(fbmTimes(z))) ^ (2 · h))
fbmPath
fbmPath(z, h) · z standard normal draws vector n · h Hurst exponent H → path vector n
Let \(z \in \mathbb{R}^n\) be standard normal draws, \(H \in (0, 1)\) (the argument h), \(C = \operatorname{fbmCovariance}(z, H)\), and \(L = \operatorname{cholesky}(C)\) the lower-triangular factor with \(C = L L^{\top}\) from Linear algebra. Then
An exact fractional Brownian path at times \(1, \dots, n\). Memory grows as \(n^2\).
Intuition Any Gaussian vector with covariance \(C\) can be made from independent normals as \(Lz\), where \(LL^{\top} = C\). The Cholesky factor mixes the independent draws so they have exactly the correlations of fractional Brownian motion, with no approximation.
Use it when you simulate long-memory processes: bursty network traffic, river levels (Hurst’s original case), rough volatility with \(H < \tfrac12\), or test data for an estimator of \(H\). Memory grows as \(n^2\) and the factorisation faster still, so stay within a few thousand points. The output is a path. Take differences to get the stationary increments (fractional Gaussian noise).
long-memorysimulation
def fbmPath(z, h) ≔ einsum("ij,j->i", cholesky(fbmCovariance(z, h)), z)
persistent
persistent(h) · h Hurst exponent H → flag 1 or 0
Let \(H \in (0, 1)\) (the argument h) be the Hurst exponent. Then
1 when increments are positively correlated.
Intuition With \(H\) above one half, increments are positively correlated: trends tend to continue and paths wander farther than Brownian motion does. Below one half, moves tend to reverse, so paths look choppy and stay closer to where they started.
Use it when you classify a fitted Hurst exponent: does the series trend (momentum, persistent river flows, long-range dependent traffic), or reverse at short scales (rough volatility)? An estimate near one half from a finite sample is not reliably either, so look at its uncertainty, not just the flag. Exactly \(H = \tfrac12\) returns 0.
long-memorytime-series
def persistent(h) ≔ h > 0.5
roughness
roughness(h) · h Hurst exponent H → fractal dimension scalar
Let \(H \in (0, 1)\) (the argument h) be the Hurst exponent. Then
The fractal dimension of the path’s graph: rougher paths have a smaller \(H\) and a larger dimension.
Intuition A path with small \(H\) is so jagged that its graph almost fills an area (dimension near 2). With \(H\) near 1 it is nearly a smooth line (dimension near 1). Brownian motion sits at 1.5.
Use it when you relate a Hurst exponent to a measured fractal dimension of a surface profile, terrain, coastline or physiological signal, or report the same property in the vocabulary a field uses. It restates \(H\) exactly, so estimate whichever is easier and convert.
long-memorymodel-interpretation
def roughness(h) ≔ 2 − h
Building blocks
fbmTimes(z)— Let \(z \in \mathbb{R}^n\). Then \(\operatorname{fbmTimes}(z)_i = i + 1\) for positions \(i = 0, \dots, n-1\), that is \(\operatorname{fbmTimes}(z) = (1, 2, \dots, n)\). The times 1, 2, …, n for a vector of length n.
Fitting an Ornstein–Uhlenbeck process
Sampled every \(\Delta t\), an Ornstein–Uhlenbeck process is an AR(1) series, \(x_{t+1} = a + b\,x_t + \varepsilon_t\) with \(b = e^{-\theta\Delta t}\). The fit is therefore one least-squares line through the pairs \((x_t, x_{t+1})\):
\[\hat\theta = -\frac{\ln \hat b}{\Delta t}, \qquad \hat\mu = \frac{\hat a}{1 - \hat b}, \qquad \hat\sigma = \frac{\operatorname{sd}(\hat\varepsilon)}{\sqrt{\Delta t}}\]The slope is floored at \(10^{-9}\) before its logarithm is taken, so a series with no positive autocorrelation reports a very large \(\hat\theta\) rather than NaN. \(\hat\sigma\) is the residual standard deviation per unit time, which matches the continuous-time \(\sigma\) when \(\theta\Delta t\) is small. The averages inside the slope and intercept run over the whole series, so on short series the estimates are biased; use a few hundred observations or more.
Throughout this section let \(x = (x_0, \dots, x_{n-1}) \in \mathbb{R}^n\) (the argument series), \(n \ge 2\), with positions counted from zero, and let
All three averages divide by \(n\), the length of the whole series, although the sums for \(\bar y\) and \(\bar x'\) have only \(n - 1\) terms; this is the source of the short-series bias. The fitted slope and intercept are \(\hat b = \operatorname{slopeFit}(x)\) and \(\hat a = \operatorname{interceptFit}(x)\), and the residuals are \(\hat\varepsilon = \operatorname{residualsOf}(x)\), all defined under the building blocks below.
thetaFrom
thetaFrom(series, dt) · series observations vector n · dt sampling interval → θ̂ scalar
Let \(x \in \mathbb{R}^n\) be the series and \(\Delta t > 0\) (the argument dt) the sampling interval. Then
The estimated reversion speed.
Intuition Sampled every \(\Delta t\), the process keeps a fraction \(b = e^{-\theta\Delta t}\) of its deviation from one observation to the next. That fraction is the slope of the line through (current value, next value), so \(\theta\) is minus the log of the slope, divided by \(\Delta t\).
Use it when you calibrate a mean-reversion speed from data: interest rates, trading spreads, commodity prices. Use a regularly sampled series of a few hundred points or more, since short series bias the slope. A slope at or below 0 is floored at \(10^{-9}\), which reports about \(20.7/\Delta t\). That means no positive autocorrelation, not fast reversion. A slope near 1 means \(\theta\) is barely identified. Pair it with muFrom and sigmaFrom, and check ouQuality.
parameter-estimationmean-reversioninterest-rate-modelling
def thetaFrom(series, dt) ≔ 0 − ln(nonzero(max(slopeFit(series), 1e-9))) / nonzero(dt)
muFrom
muFrom(series) · series observations vector n → μ̂ scalar
Let \(x \in \mathbb{R}^n\) be the series. Then
\[\operatorname{muFrom}(x) = \hat\mu = \frac{\hat a}{\operatorname{nonzero}(1 - \hat b)} .\]The estimated long-run mean.
Intuition The fitted line \(x_{t+1} = a + b\,x_t\) has one fixed point, where the next value equals the current one: \(a/(1-b)\). That is the level the series is pulled toward.
Use it when you estimate the long-run level of a mean-reverting series alongside thetaFrom and sigmaFrom. It is fragile when the slope is close to 1 (slow reversion): \(1 - \hat b\) is then small, and small errors in the intercept swing the estimate a lot. In that case treat it with caution, or compare it with the plain average of a long series. Like the rest of this fit, it is biased on short series.
parameter-estimationmean-reversion
def muFrom(series) ≔ interceptFit(series) / nonzero(1 − slopeFit(series))
sigmaFrom
sigmaFrom(series, dt) · series observations vector n · dt sampling interval → σ̂ scalar
Let \(x \in \mathbb{R}^n\) be the series, \(\Delta t > 0\) (the argument dt), \(\hat\varepsilon = \operatorname{residualsOf}(x) \in \mathbb{R}^n\) and \(\bar\varepsilon = \frac{1}{n} \sum_{i=0}^{n-1} \hat\varepsilon_i\). Then
The estimated volatility, the sample standard deviation of the one-step residuals over \(\sqrt{\Delta t}\). The sum includes the zero residual at the last position.
Intuition Whatever the line fails to predict one step ahead is the noise that came in during that step. Its size grows like \(\sigma\sqrt{\Delta t}\), so dividing by \(\sqrt{\Delta t}\) gives the volatility per unit time.
Use it when you complete a calibration from a series. It matches the continuous-time \(\sigma\) only when \(\theta\,\Delta t\) is small. For coarse sampling, the exact relation is \(\sigma = \operatorname{sd} \cdot \sqrt{2\theta / (1 - e^{-2\theta\Delta t})}\), which you can apply with thetaFrom. On short series it shares the bias of the rest of the fit.
parameter-estimationmean-reversion
def sigmaFrom(series, dt) ≔ std(residualsOf(series)) / nonzero(sqrt(dt))
ouQuality
ouQuality(series) · series observations vector n → R² scalar
Let \(x \in \mathbb{R}^n\) be the series, \(\hat\varepsilon = \operatorname{residualsOf}(x)\), and \(y = \operatorname{laterOf}(x) = (x_1, \dots, x_{n-1}, 0)\). Write \(s^2(v) = \frac{1}{n-1} \sum_{i=0}^{n-1} (v_i - \bar v)^2\) for the sample variance of a vector \(v \in \mathbb{R}^n\). Then
\[\operatorname{ouQuality}(x) = 1 - \frac{s^2(\hat\varepsilon)}{\operatorname{nonzero}\big(s^2(y)\big)} .\]The share of the variance of the next value that the line explains, \(1 - \operatorname{Var}(\hat\varepsilon)/\operatorname{Var}(x_{t+1})\). Both variances run over all \(n\) positions, including the zero placed in the last one; for a series whose level is far from zero that zero inflates \(s^2(y)\), so the reported value is higher than the \(R^2\) of an exact regression.
Intuition The share of the next value’s variance that the current value explains. A high value means the series is very predictable one step ahead, so reversion is slow compared with the sampling interval. A low value means noise dominates.
Use it when you sanity-check a fit before trusting \(\hat\theta\) and \(\hat\mu\). A low value together with a slope near 0 means there is little reversion signal to estimate. It does not test whether an Ornstein–Uhlenbeck model is the right one (inspect the residuals for that). For series far from zero, the zero in the last position inflates it, so use it to compare fits of similar series, not as an exact \(R^2\).
parameter-estimationmodel-validation
def ouQuality(series) ≔ 1 − variance(residualsOf(series)) / nonzero(variance(laterOf(series) · paired(series)))
Building blocks
laterOf(series)— Let \(x = (x_0, \dots, x_{n-1}) \in \mathbb{R}^n\). Then \(\operatorname{laterOf}(x)_i = \sum_{j=0}^{n-1} [\, i - j = -1 \,]\, x_j\), which is \(x_{i+1}\) for \(i \le n-2\) and \(0\) for \(i = n-1\). The series shifted forward one step: entry t holds the value at t + 1, and the last entry is 0.paired(series)— Let \(x \in \mathbb{R}^n\). Then \(\operatorname{paired}(x)_i = [\, i < n - 1 \,]\) for \(i = 0, \dots, n-1\). 1 for every position that has a next value, 0 for the last.slopeFit(series)— Let \(x \in \mathbb{R}^n\) with \(\bar x\) and \(\bar y\) as above. Then \(\operatorname{slopeFit}(x) = \hat b = \sum_{i=0}^{n-2} (x_i - \bar x)(x_{i+1} - \bar y) \,\big/\, \operatorname{nonzero}\big( \sum_{i=0}^{n-2} (x_i - \bar x)^2 \big)\). The least-squares slope b̂ of the next value on the current one.interceptFit(series)— Let \(x \in \mathbb{R}^n\) with \(\bar y\), \(\bar x'\) as above and \(\hat b = \operatorname{slopeFit}(x)\). Then \(\operatorname{interceptFit}(x) = \hat a = \bar y - \hat b\, \bar x'\). The least-squares intercept â.residualsOf(series)— Let \(x \in \mathbb{R}^n\), \(\hat a = \operatorname{interceptFit}(x)\) and \(\hat b = \operatorname{slopeFit}(x)\). Then \(\operatorname{residualsOf}(x)_i = x_{i+1} - \hat a - \hat b\, x_i\) for \(i \le n-2\), and \(\operatorname{residualsOf}(x)_{n-1} = 0\). The one-step residuals of the fitted line, zero at the last position.
Examples
One Euler and one Milstein step
State 1, drift 0.5, diffusion 0.2 with derivative 0.2 (as for GBM with \(\sigma = 0.2\)), time step 0.01 and Brownian increment 0.2.
open Prelude.Diffusions
fn eulerStep(x, a, b, dt, dw) ↦ euler(x, a, b, dt, dw)
fn milsteinStep(x, a, b, db, dt, dw) ↦ milstein(x, a, b, db, dt, dw)
using Examples;
var e = Sde.EulerStep(1.0, 0.5, 0.2, 0.01, 0.2);
var m = Sde.MilsteinStep(1.0, 0.5, 0.2, 0.2, 0.01, 0.2);
Output (verified):
euler = 1.045
milstein = 1.0456
Euler adds \(0.5 \cdot 0.01 + 0.2 \cdot 0.2 = 0.045\); Milstein adds \(\tfrac12 \cdot 0.2 \cdot 0.2 \cdot (0.04 - 0.01) = 0.0006\) more.
The mean and median of a stock price
A price of 100 with growth 5% and volatility 20% over one year.
open Prelude.Diffusions
fn expectedPrice(x0, mu, t) ↦ gbmMean(x0, mu, t)
fn typicalPrice(x0, mu, sigma, t) ↦ gbmMedian(x0, mu, sigma, t)
fn exactPrice(x0, mu, sigma, t, z) ↦ gbmExact(x0, mu, sigma, t, z)
using Examples;
var mean = Sde.ExpectedPrice(100.0, 0.05, 1.0);
var median = Sde.TypicalPrice(100.0, 0.05, 0.2, 1.0);
var upOneSd = Sde.ExactPrice(100.0, 0.05, 0.2, 1.0, 1.0);
Output (verified, rounded to 6 decimal places):
gbmMean = 105.12711
gbmMedian = 103.045453
gbmExact = 125.860001
The mean is \(100e^{0.05}\) but the median only \(100e^{0.03}\). A draw one standard deviation up gives \(100e^{0.03 + 0.2}\).
Mean reversion and the Feller condition
open Prelude.Diffusions
fn settledVariance(theta, sigma) ↦ ouStationary(theta, sigma)
fn halfway(theta) ↦ halfLife(theta)
fn revertedMean(theta, mu, x0, t) ↦ ouMean(theta, mu, x0, t)
fn positiveRate(theta, mu, sigma) ↦ fellerHolds(theta, mu, sigma)
using Examples;
var variance = Sde.SettledVariance(2.0, 1.0);
var half = Sde.Halfway(2.0);
var meanAtHalf = Sde.RevertedMean(2.0, 1.0, 0.0, 0.5);
var fastRates = Sde.PositiveRate(2.0, 0.04, 0.3);
var slowRates = Sde.PositiveRate(1.0, 0.04, 0.3);
Output (verified, rounded to 6 decimal places):
ouStationary = 0.25
halfLife = 0.346574
ouMean = 0.632121
feller = 1
feller2 = 0
With \(\theta = 2\) and \(\sigma = 1\) the stationary variance is \(1/4\) and the half-life \(\ln 2/2\). Starting at 0 with mean 1, the expected value at \(t = 0.5\) is \(1 - e^{-1}\). A rate model with mean 4% and volatility 0.3 satisfies Feller’s condition at \(\theta = 2\) (\(0.16 \ge 0.09\)) but not at \(\theta = 1\) (\(0.08 < 0.09\)).
Quadratic variation and fractional covariance
open Prelude.Diffusions
fn squaredMoves(path) ↦ quadraticVariation(path)
fn fbmCov(z, h) ↦ fbmCovariance(z, h)
fn fatTail(sigma, rate, jumpSpread, t) ↦ mertonKurtosis(sigma, rate, jumpSpread, t)
using Examples;
using MLambda.Hilbert.Runtime;
var qv = Sde.SquaredMoves((Tensor)new[] { 0d, 1, 0, 2 });
var cov = Sde.FbmCov((Tensor)new[] { 0d, 0, 0 }, 0.5);
var kurt = Sde.FatTail(0.2, 1.0, 0.1, 1.0);
Output (verified):
qv = 6
fbmCov = shape [3, 3] values [1, 1, 1, 1, 2, 2, 1, 2, 3]
kurtosis = 3.12
The increments are 0, 1, −1 and 2, whose squares sum to 6. At \(H = \tfrac12\) the covariance reduces to \(\min(s, t)\), the covariance of ordinary Brownian motion. One jump a year of spread 0.1 on top of 20% volatility adds \(3 \cdot 10^{-4} / 0.05^2 = 0.12\) to the kurtosis.
Recovering Ornstein–Uhlenbeck parameters from a series
The series is 400 steps of \(x_{t+1} = 1 + 0.5\,x_t + 0.1\,\varepsilon_t\) with \(\Delta t = 1\), so the true values are \(\theta = \ln 2 \approx 0.693\), \(\mu = 2\) and a one-step noise of 0.1.
open Prelude.Diffusions
fn fittedTheta(series, dt) ↦ thetaFrom(series, dt)
fn fittedMu(series) ↦ muFrom(series)
fn fittedSigma(series, dt) ↦ sigmaFrom(series, dt)
fn fittedQuality(series) ↦ ouQuality(series)
using Examples;
using MLambda.Hilbert.Runtime;
var rng = new System.Random(7);
var x = new double[400];
x[0] = 2;
for (int i = 1; i < x.Length; i++)
{
double g = Math.Sqrt(-2 * Math.Log(1 - rng.NextDouble())) * Math.Cos(2 * Math.PI * rng.NextDouble());
x[i] = 1 + 0.5 * x[i - 1] + 0.1 * g;
}
var theta = Sde.FittedTheta((Tensor)x, 1.0);
var mu = Sde.FittedMu((Tensor)x);
var sigma = Sde.FittedSigma((Tensor)x, 1.0);
var quality = Sde.FittedQuality((Tensor)x);
Output (verified, rounded to 6 decimal places):
theta = 0.725596
mu = 1.987928
sigma = 0.098042
quality = 0.572377
An exact least-squares regression on the same 400 pairs gives \(\hat\theta = 0.7258\) and \(\hat\mu = 1.9929\), so the difference from the true values is sampling error rather than the estimator. On a six-point series the same functions are noticeably biased, as the section above explains.