A step's volatility: the spread grows as the root of the time¶
A random walk's standard deviation over a step, from its standard deviation over a period, sigma_dt = sigma_T sqrt(dt / T).
Independent increments add their variances, so the variance grows as the time and the spread as its root: n steps of dt = T / n give back the period's variance, sigma_T^2. Both spreads are the walk's own, dimensionless (a log price's, the spread of a ratio); the period is the game's (Bourse's minute, 60 s), never a literal, so nothing is a volatility per root second. Source: Hull, Options, Futures and Other Derivatives, 15.1 (the volatility over a time scales as its root). Checked against a reference in numpy, which has its own tests.
- Source
- Hull, Options, Futures and Other Derivatives, 15.1.
- Proof
- curve, invariant, limit.
step_volatility.trait
namespace lusine.stochastic;
float step_volatility(float per_period, float<s> period, float<s> dt)
{
return per_period * sqrt(dt / period);
}
A random walk's standard deviation over a step dt, from its standard deviation over a period T, σ_dt = σ_T √(dt / T). Independent increments add their variances, so n steps of T / n give back the period's variance, σ_T². Both spreads are dimensionless (a log price's, as Bourse states its factors' a minute); the period is the game's, passed in.
import os, sys
root = os.getcwd()
while not os.path.isfile(os.path.join(root, 'docsite.json')):
root = os.path.dirname(root)
sys.path.insert(0, os.path.join(root, 'tools'))
import numpy as np
import matplotlib.pyplot as plt
from lusine import trait_run
EPS = float(np.finfo(np.float32).eps) # Trait's float is binary32
U = EPS / 2 # binary32's unit roundoff, 2^-24
f32 = lambda a: np.float32(a).astype(float)
The kernel¶
The law as the library states it.
print(trait_run.source('stochastic/step_volatility.trait'))
// lusine.stochastic.step_volatility: a random walk's standard deviation over a step, from its standard deviation over a
// period, sigma_dt = sigma_T sqrt(dt / T). Independent increments add their variances, so the variance grows as the
// time and the spread as its root: n steps of dt = T / n give back the period's variance, sigma_T^2. Both spreads are
// the walk's own, dimensionless (a log price's, the spread of a ratio); the period is the
// game's (Bourse's minute, 60 s), never a literal, so nothing is a volatility per root second.
// Source: Hull, Options, Futures and Other Derivatives, 15.1 (the volatility over a time scales as its root).
// Checked against a reference in numpy, which has its own tests.
namespace lusine.stochastic;
float step_volatility(float per_period, float<s> period, float<s> dt)
{
return per_period * sqrt(dt / period);
}
The program¶
Each section runs the kernel over a lattice of items, one a case: the host sets the spread, the period and the step, and main() binds the kernel over the items.
print(trait_run.law_program('stochastic/step_volatility.trait', 'step_volatility'))
def vol(per_period, period, dt):
return trait_run.law('stochastic/step_volatility.trait', 'step_volatility', per_period=per_period, period=period,
dt=dt)
namespace probe;
lusine.version version = { { number = 1 } };
struct item {
host float per_period;
host float<s> period;
host float<s> dt;
float result_0;
};
lattice<item, 1> items(spacing = 1 m, boundary = wall);
void main() { items.result_0 = lusine.stochastic.step_volatility(items.per_period, items.period, items.dt); }
Curve: σ = 0.01 a minute is 0.0775 an hour¶
The textbook's case: a spread of 0.01 a minute is 0.01 √60 ≈ 0.0775 over an hour, and 0.01 / √60 over a second. On a log-log plot the law is a line of slope ½; the planted fault, the spread scaled linearly (dt / T without the root), is a line of slope 1, and meets it only at dt = T.
v = vol(0.01, 60.0, [3600.0, 1.0])
assert abs(v[0] - 0.01 * np.sqrt(60.0)) <= 2 * EPS * v[0] and round(v[0], 4) == 0.0775, v
assert abs(v[1] - 0.01 / np.sqrt(60.0)) <= 2 * EPS * v[1], v
print(f'an hour: {v[0]:.6f}; a second: {v[1]:.6g}')
dt = f32(np.geomspace(0.1, 3600.0, 60))
v = vol(0.01, 60.0, dt)
slope = np.polyfit(np.log(dt), np.log(v), 1)[0]
assert abs(slope - 0.5) < 1e-6, slope
fig, ax = plt.subplots(figsize=(7, 3.2))
ax.loglog(dt, v, label='step_volatility, σ_T √(dt / T)')
ax.loglog(dt, 0.01 * dt / 60.0, '--', label='planted fault: σ_T dt / T')
ax.axvline(60.0, color='grey', lw=0.5)
ax.set_xlabel('dt (s)')
ax.set_ylabel('σ over the step')
ax.legend()
plt.show()
an hour: 0.077460; a second: 0.00129099
Invariant: n steps' variance is the period's¶
Over n steps of dt = T / n the walk's variance is n σ_dt², which must be σ_T² at any n. At n a power of four the root is a power of two and the step's spread exact in binary32, so the variance comes back bit for bit; at any other n it comes back within 6 × 2⁻²⁴, the square doubling the spread's own 2.5 × 2⁻²⁴ (the last section). The same holds over a session: a factor's variance over 6.5 hours, (S / dt) σ_dt², is the same at any dt, as Bourse's clock asks.
sigma, T = f32(0.012), 60.0
n = 4.0 ** np.arange(0, 9)
v = vol(sigma, T, T / n)
assert np.array_equal(n * v * v, np.full(n.shape, sigma * sigma)), n * v * v
n = np.arange(1, 1001, dtype=float)
dt = f32(T / n)
v = vol(sigma, T, dt)
var = (T / dt) * v * v
assert np.all(np.abs(var / (sigma * sigma) - 1) <= 6 * U), np.max(np.abs(var / (sigma * sigma) - 1)) / U
print(f'n sigma_dt^2 / sigma_T^2 - 1 over n = 1..1000: at most {np.max(np.abs(var / (sigma * sigma) - 1)) / U:.2f} x 2^-24')
S = 23400.0
dt = f32([0.25, 1.0, 15.0, 60.0, 300.0, 1800.0])
session = S / dt * vol(sigma, T, dt) ** 2
want = S / T * sigma * sigma
assert np.all(np.abs(session / want - 1) <= 6 * U), session
print('a session\'s variance at dt =', dt, ':', session, 'want', want)
n sigma_dt^2 / sigma_T^2 - 1 over n = 1..1000: at most 3.68 x 2^-24 a session's variance at dt = [2.5e-01 1.0e+00 1.5e+01 6.0e+01 3.0e+02 1.8e+03] : [0.05616 0.05616 0.05616 0.05616 0.05616001 0.05616001] want 0.056160000976324086
Limit: a step of the period, and of nothing¶
A step of the period is the period's spread, bit for bit (dt / T is 1, its root 1); a step of nothing has no spread. A game refuses a step of nothing before it gets here (Bourse's clock does), but the kernel's answer at the limit is the right one.
s = f32([0.0, 1e-6, 0.01, 0.35, 2.0])
assert np.array_equal(vol(s, 60.0, 60.0), s) and np.array_equal(vol(s, 7.5, 7.5), s)
assert np.all(vol(s, 60.0, 0.0) == 0.0)
print('dt = T gives sigma_T itself; dt = 0 gives 0')
dt = T gives sigma_T itself; dt = 0 gives 0
Agreement: binary32 against the float64 reference¶
Three thousand cases, spreads from 10⁻⁴ to 1 and steps from a sixtieth of a second to four periods, each input a binary32: the emitted kernel against the reference in float64 on the same inputs. A quotient, a correctly rounded root and a product bound the relative error by 2.5 × 2⁻²⁴; that is what it is held to.
from lusine.kernels.stochastic.step_volatility import step_volatility as reference
rng = np.random.default_rng(302)
s = f32(10 ** rng.uniform(-4, 0, 3000))
dt = f32(60 * 10 ** rng.uniform(-3.5, 0.6, 3000))
rel = np.abs(vol(s, 60.0, dt) / reference(s, 60.0, dt) - 1)
assert np.max(rel) <= 2.5 * U, np.max(rel) / U
print(f'relative error at most {np.max(rel):.3g} = {np.max(rel) / U:.2f} x 2^-24')
relative error at most 1.19e-07 = 2.00 x 2^-24