diff --git a/docs/theory/levy.md b/docs/theory/levy.md index fdf42710..702c0663 100644 --- a/docs/theory/levy.md +++ b/docs/theory/levy.md @@ -41,7 +41,9 @@ While $\tau_t$ is always continuous, $\lambda$ can exhibit jumps. Since the time \Phi_{y_t, u} = {\mathbb E}\left[e^{i u x_{\tau_t}}\right] = {\mathbb E}\left[{\mathbb E}\left[\left.e^{i u x_s}\right|\tau_t=s\right]\right] \end{equation} -where the inside expectation is taken on $x_{\tau_t}$ conditional on a fixed value of $\tau_t = s$ and the outside expectation is on all possible values of $\tau_t$. If the random time $\tau_t$ is independent of $x_t$, the randomness due to the Lévy process can be integrated out using the characteristic function of $x_t$: +where the inside expectation is taken on $x_{\tau_t}$ conditional on a fixed value of $\tau_t = s$ and the outside expectation is on all possible values of $\tau_t$. + +If the random time $\tau_t$ is independent of $x_t$, the randomness due to the Lévy process can be integrated out using the characteristic function of $x_t$: \begin{equation} \Phi_{y_t, u} = {\mathbb E}\left[e^{-\tau_t \phi_{x,u}}\right] = {\mathbb L}_{\tau_t}\left(\phi_{x_1,u}\right) diff --git a/quantflow/sp/base.py b/quantflow/sp/base.py index 896087ec..03c6ef9c 100755 --- a/quantflow/sp/base.py +++ b/quantflow/sp/base.py @@ -204,7 +204,7 @@ def integrated_log_laplace( r"""The log-Laplace transform of the cumulative process: \begin{equation} - e^{\phi_{t, u}} = {\mathbb E} \left[e^{i u \int_0^t x_s ds}\right] + \phi_{t, u} = \log {\mathbb E} \left[e^{-u \int_0^t x_s ds}\right] \end{equation} """ diff --git a/quantflow/sp/bns.py b/quantflow/sp/bns.py index a7fa1f35..2cf66674 100644 --- a/quantflow/sp/bns.py +++ b/quantflow/sp/bns.py @@ -101,18 +101,29 @@ def sample(self, n: int, time_horizon: float = 1, time_steps: int = 100) -> Path return self.sample_from_draws(Paths.normal_draws(n, time_horizon, time_steps)) def sample_from_draws(self, path_dw: Paths, *args: Paths) -> Paths: + kappa = self.variance_process.kappa + time_steps = path_dw.time_steps + time_horizon = path_dw.t + dt = path_dw.dt + if args: path_dz = args[0] else: - # generate the background driving process samples if not provided + # BDLP runs at kappa-rescaled time, so sample over [0, kappa*T] path_dz = self.variance_process.bdlp.sample( - path_dw.samples, path_dw.t, path_dw.time_steps + path_dw.samples, kappa * time_horizon, time_steps ) - dt = path_dw.dt - # sample the activity rate process - v = self.variance_process.sample_from_draws(path_dz) - # create the time-changed Brownian motion - dw = path_dw.data * np.sqrt(v.data * dt) - paths = np.zeros(dw.shape) - paths[1:] = np.cumsum(dw[:-1], axis=0) + path_dz.data - return Paths(t=path_dw.t, data=paths) + + # Variance via the OU recursion v_{i+1} = v_i * e^{-kappa dt} + Delta z_{i+1} + decay = np.exp(-kappa * dt) + dz = np.diff(path_dz.data, axis=0) + v = np.zeros_like(path_dz.data) + v[0] = self.variance_process.rate + for i in range(time_steps): + v[i + 1] = v[i] * decay + dz[i] + + # Price: x_t = integral sqrt(v) dW + rho * z_{kappa t} + diffusion = np.sqrt(v[:-1] * dt) * path_dw.data[:-1] + paths = np.zeros_like(path_dw.data) + paths[1:] = np.cumsum(diffusion, axis=0) + self.rho * path_dz.data[1:] + return Paths(t=time_horizon, data=paths) diff --git a/quantflow/sp/dsp.py b/quantflow/sp/dsp.py index c18b404b..0c4a98ed 100755 --- a/quantflow/sp/dsp.py +++ b/quantflow/sp/dsp.py @@ -34,8 +34,10 @@ def support(self, mean: float, std: float, points: int) -> FloatArray: return np.linspace(0, points, points + 1) def characteristic_exponent(self, t: FloatArrayLike, u: Vector) -> Vector: - phi = self.poisson.characteristic_exponent(t, u) - return -self.intensity.integrated_log_laplace(t, phi) + # phi_x is the per-unit-time Lévy exponent of the inner Poisson; + # the integrated log-Laplace already absorbs the time horizon + phi_x = self.poisson.characteristic_exponent(1, u) + return -self.intensity.integrated_log_laplace(t, phi_x) def arrivals(self, t: float = 1) -> FloatArray: paths = self.intensity.sample(1, t, math.ceil(100 * t)).integrate() diff --git a/quantflow/sp/heston.py b/quantflow/sp/heston.py index 1131fa1c..f1c83e17 100644 --- a/quantflow/sp/heston.py +++ b/quantflow/sp/heston.py @@ -315,6 +315,13 @@ def characteristic_exponent( t, u ) + self.jumps.characteristic_exponent(t, u) + def sample_from_draws(self, path1: Paths, *args: Paths) -> Paths: + diffusion = super().sample_from_draws(path1, *args) + jump_path = self.jumps.sample( + diffusion.samples, diffusion.t, diffusion.time_steps + ) + return Paths(t=diffusion.t, data=diffusion.data + jump_path.data) + class DoubleHeston(StochasticProcess1D): r"""Double Heston stochastic volatility model. diff --git a/quantflow/sp/ou.py b/quantflow/sp/ou.py index 908f0068..dde0fe97 100755 --- a/quantflow/sp/ou.py +++ b/quantflow/sp/ou.py @@ -5,24 +5,25 @@ import numpy as np from pydantic import Field -from scipy.optimize import Bounds from scipy.stats import gamma, norm from ..ta.paths import Paths from ..utils.distributions import Exponential from ..utils.types import Float, FloatArrayLike, Vector -from .base import Im, IntensityProcess +from .base import IntensityProcess, StochasticProcess1D from .poisson import CompoundPoissonProcess, D from .weiner import WeinerProcess -class Vasicek(IntensityProcess): - r"""Gaussian OU process, also know as the - [Vasiceck model](https://en.wikipedia.org/wiki/Vasicek_model) +class Vasicek(StochasticProcess1D): + r"""Gaussian OU process, also known as the + [Vasicek model](https://en.wikipedia.org/wiki/Vasicek_model). Historically, the Vasicek model was used to model the short rate, but it can be used to model any process that reverts to a mean level at a rate proportional to - the difference between the current level and the mean level. + the difference between the current level and the mean level. Unlike intensity + processes, the Vasicek process can take negative values, so it is not constrained + to a positive domain. \begin{equation} dx_t = \kappa (\theta - x_t) dt + \sigma dw_t @@ -37,11 +38,15 @@ class Vasicek(IntensityProcess): \end{equation} """ + rate: float = Field(default=1.0, description=r"Initial value $x_0$") + kappa: float = Field( + default=1.0, gt=0, description=r"Mean reversion speed $\kappa$" + ) + theta: float = Field(default=1.0, description=r"Mean level $\theta$") bdlp: WeinerProcess = Field( default_factory=WeinerProcess, - description="Background driving Weiner process", + description="Background driving Wiener process", ) - theta: float = Field(default=1.0, gt=0, description=r"Mean rate $\theta$") def characteristic_exponent(self, t: FloatArrayLike, u: Vector) -> Vector: r"""The characteristic exponent of the Vasicek process is given by @@ -54,9 +59,6 @@ def characteristic_exponent(self, t: FloatArrayLike, u: Vector) -> Vector: var = self.analytical_variance(t) return u * (-1j * mu + 0.5 * var * u) - def integrated_log_laplace(self, t: FloatArrayLike, u: Vector) -> Vector: - raise NotImplementedError - def sample(self, n: int, time_horizon: float = 1, time_steps: int = 100) -> Paths: paths = Paths.normal_draws(n, time_horizon, time_steps) return self.sample_from_draws(paths) @@ -74,8 +76,8 @@ def sample_from_draws(self, draws: Paths, *args: Paths) -> Paths: paths[t + 1, :] = x + dx return Paths(t=draws.t, data=paths) - def domain_range(self) -> Bounds: - return Bounds(-np.inf, np.inf) + def ekt(self, t: FloatArrayLike) -> FloatArrayLike: + return np.exp(-self.kappa * t) def analytical_mean(self, t: FloatArrayLike) -> FloatArrayLike: r"""The analytical mean of the Vasicek process is given by @@ -149,10 +151,12 @@ def intensity(self) -> float: @abstractmethod def characteristic_exponent(self, t: FloatArrayLike, u: Vector) -> Vector: - r""" + r"""Characteristic exponent of a non-Gaussian OU process driven by a + Lévy process $Z$ with characteristic exponent $\phi_Z$: + \begin{equation} - \phi_{x_t,u} = iu x_0 e^{-\kappa t} + - \int_{ue^{-\kappa t}}^{u} \frac{\psi_Z(v)}{v} , dv + \phi_{x_t, u} = - iu\, x_0\, e^{-\kappa t} + + \int_{u e^{-\kappa t}}^{u} \frac{\phi_Z(v)}{v}\,dv \end{equation} """ @@ -180,30 +184,82 @@ def create(cls, rate: float = 1, decay: float = 1, kappa: float = 1) -> Self: ) def characteristic_exponent(self, t: FloatArrayLike, u: Vector) -> Vector: - r""" + r"""Closed form of the [NGOU][quantflow.sp.ou.NGOU] characteristic + exponent when the BDLP is a compound Poisson process with intensity + $\lambda$ and exponential jumps of decay $\beta$: + + \begin{equation} + \phi_{x_t, u} = + - iu\, x_0\, e^{-\kappa t} + + \lambda \log\!\left( + \frac{\beta - iu}{\beta - iu\, e^{-\kappa t}} + \right) + \end{equation} + + Derivation. The characteristic function of the + [Exponential][quantflow.utils.distributions.Exponential.characteristic] + jumps is + \begin{equation} - \phi_{x_t,u} = - - iu x_0 e^{-\kappa t} - + \lambda \log\!\left(\frac{\beta - iu}{\beta - iu e^{-\kappa t}}\right) + \Phi_J(v) = \frac{\beta}{\beta - iv}, \end{equation} + + so the BDLP unit-time characteristic exponent + ([CompoundPoissonProcess][quantflow.sp.poisson.CompoundPoissonProcess] + at $t = 1$) is + + \begin{equation} + \phi_Z(v) = \lambda\bigl(1 - \Phi_J(v)\bigr) + = -\frac{i v\, \lambda}{\beta - iv}, + \end{equation} + + and $\phi_Z(v)/v = -i\lambda / (\beta - iv)$. Substituting + $w = \beta - iv$ (so $dv = i\, dw$) in the NGOU integral gives + + \begin{equation} + \int_{u e^{-\kappa t}}^{u} \frac{\phi_Z(v)}{v}\, dv + = \lambda \int_{\beta - iu\, e^{-\kappa t}}^{\beta - iu} + \frac{dw}{w} + = \lambda \log\!\left( + \frac{\beta - iu}{\beta - iu\, e^{-\kappa t}} + \right), + \end{equation} + + which combined with the drift term $-iu\, x_0\, e^{-\kappa t}$ yields + the closed form above. """ b = self.beta - iu = Im * u + iu = 1j * u c1 = iu * np.exp(-self.kappa * t) c0 = self.intensity * np.log((b - c1) / (b - iu)) return -c0 - c1 * self.rate def integrated_log_laplace(self, t: FloatArrayLike, u: Vector) -> Vector: + r"""Log-Laplace transform of the time-integrated Gamma-OU process. + + \begin{equation} + \begin{aligned} + \phi(t, u) &= \log E\!\left[e^{-u \int_0^t x_s\, ds}\right] \\ + &= -u\, x_0\, \frac{1 - e^{-\kappa t}}{\kappa} + + \frac{\lambda}{u + \beta\kappa} + \!\left[\beta\kappa\, + \log\!\frac{\beta\kappa + u(1 - e^{-\kappa t})}{\beta\kappa} + - u\kappa t\right] + \end{aligned} + \end{equation} + + where $x_0$ is the initial rate, $\lambda$ is the BDLP intensity and + $\beta$ the exponential jump decay. + """ kappa = self.kappa b = self.beta - iu = Im * u - iuk = iu / kappa + muk = -u / kappa ekt = np.exp(-kappa * t) - c1 = iuk * (1 - ekt) + c1 = muk * (1 - ekt) c0 = self.intensity * ( - b * np.log(b / (iuk + (b - iuk) / ekt)) / (iuk - b) - kappa * t + b * np.log(b / (muk + (b - muk) / ekt)) / (muk - b) - kappa * t ) - return np.exp(c0 + c1 * self.rate) + return c0 + c1 * self.rate def sample(self, n: int, time_horizon: float = 1, time_steps: int = 100) -> Paths: dt = time_horizon / time_steps @@ -235,28 +291,43 @@ def _advance( ) -> int: x = pp[i - 1] kappa = self.kappa - t0 = i * dt - t1 = t0 + dt + # step i fills pp[i] from pp[i-1] over the interval [(i-1) dt, i dt] + t0 = (i - 1) * dt + t1 = i * dt a = arrival or t1 pp[i] = x - kappa * x * (a - t0) - kappa * (x + jump) * (t1 - a) + jump return i + 1 - def cumulative_characteristic2(self, t: FloatArrayLike, u: Vector) -> Vector: - """Formula from a paper""" - kappa = self.kappa - b = self.beta - iu = Im * u - iuk = iu / kappa - ekt = np.exp(-kappa * t) - c1 = iuk * (1 - ekt) - c0 = self.intensity * (b * np.log(b / (b - c1)) - iu * t) / (iuk - b) - return np.exp(c0 + c1 * self.rate) - def analytical_mean(self, t: FloatArrayLike) -> FloatArrayLike: - return self.intensity / self.beta + r"""Analytical mean of the Gamma-OU process at time $t$. + + \begin{equation} + E[x_t] = x_0\, e^{-\kappa t} + + \frac{\lambda}{\beta}\bigl(1 - e^{-\kappa t}\bigr) + \end{equation} + + which converges to the stationary mean $\lambda / \beta$ as + $t \to \infty$. + """ + ekt = self.ekt(t) + return self.rate * ekt + (self.intensity / self.beta) * (1 - ekt) def analytical_variance(self, t: FloatArrayLike) -> FloatArrayLike: - return self.intensity / self.beta / self.beta + r"""Analytical variance of the Gamma-OU process at time $t$. + + \begin{equation} + \mathrm{Var}(x_t) = \frac{\lambda}{\beta^2} + \bigl(1 - e^{-2\kappa t}\bigr) + \end{equation} + + which converges to the stationary variance $\lambda / \beta^2$ as + $t \to \infty$. + """ + return self.intensity * (1 - self.ekt(2 * t)) / (self.beta * self.beta) def analytical_pdf(self, t: FloatArrayLike, x: FloatArrayLike) -> FloatArrayLike: + """Stationary marginal density of the Gamma-OU process. The transient + density is not available in closed form; use `pdf_from_characteristic` + for finite $t$. + """ return gamma.pdf(x, self.intensity, scale=1 / self.beta) diff --git a/quantflow_tests/test_bns.py b/quantflow_tests/test_bns.py index 92cf4f10..096d1938 100644 --- a/quantflow_tests/test_bns.py +++ b/quantflow_tests/test_bns.py @@ -32,10 +32,10 @@ def test_bns_horizon(bns: BNS) -> None: def test_bns_no_leverage_matches_integrated_laplace(bns: BNS) -> None: # at rho=0: E[exp(iu x_t)] = E[exp(-u^2/2 * int v_s ds)] - # which is the GammaOU integrated Laplace transform at i*u^2/2 - u = np.array([0.1, 0.5, 1.0, 2.0, 5.0], dtype=complex) + # which is the GammaOU integrated Laplace transform evaluated at u^2/2 + u = np.array([0.1, 0.5, 1.0, 2.0, 5.0]) cf_bns = bns.characteristic(1, u) - cf_via_ou = bns.variance_process.integrated_log_laplace(1, 1j * u * u / 2) + cf_via_ou = np.exp(bns.variance_process.integrated_log_laplace(1, u * u / 2)) np.testing.assert_allclose(cf_bns, cf_via_ou, atol=1e-12) @@ -45,3 +45,16 @@ def test_bns_leverage(bns_leverage: BNS) -> None: assert bns_leverage.characteristic(1, 0) == 1 # E[x_t] = rho * kappa * t * intensity / beta = -0.5 * 1 * 1 * 1.25 / 5 assert m.mean_from_characteristic() == pytest.approx(-0.125, 1e-3) + + +def test_bns_sample_moments(bns_leverage: BNS) -> None: + np.random.seed(42) + paths = bns_leverage.sample(5000, time_horizon=1.0, time_steps=200) + assert paths.data.shape == (201, 5000) + # Compare empirical moments at T=1 against the characteristic-function moments + m = bns_leverage.marginal(1) + terminal = paths.data[-1] + assert terminal.mean() == pytest.approx( + float(m.mean_from_characteristic()), abs=0.02 + ) + assert terminal.std() == pytest.approx(float(m.std_from_characteristic()), abs=0.02) diff --git a/quantflow_tests/test_heston.py b/quantflow_tests/test_heston.py index de36bc12..71f68722 100644 --- a/quantflow_tests/test_heston.py +++ b/quantflow_tests/test_heston.py @@ -1,3 +1,4 @@ +import numpy as np import pytest from quantflow.sp.heston import DoubleHeston, DoubleHestonJ, Heston, HestonJ @@ -83,3 +84,21 @@ def test_double_heston_jumps_characteristic( characteristic_tests(m) assert m.mean() == pytest.approx(0.0, abs=1e-6) assert m.std() == pytest.approx(0.5, rel=0.05) + + +def test_heston_jumps_sample_includes_jumps( + heston_jumps: HestonJ[DoubleExponential], +) -> None: + np.random.seed(42) + paths = heston_jumps.sample(5000, time_horizon=1.0, time_steps=200) + # If jumps were dropped, std would collapse to sqrt(1 - jump_fraction) * 0.5 + # ~ 0.418, well below 0.5 + assert paths.data[-1].std() == pytest.approx(0.5, abs=0.03) + + +def test_double_heston_jumps_sample_includes_jumps( + double_heston_jumps: DoubleHestonJ[DoubleExponential], +) -> None: + np.random.seed(42) + paths = double_heston_jumps.sample(5000, time_horizon=1.0, time_steps=200) + assert paths.data[-1].std() == pytest.approx(0.5, abs=0.03) diff --git a/quantflow_tests/test_ou.py b/quantflow_tests/test_ou.py index d3ef0759..62a822fe 100644 --- a/quantflow_tests/test_ou.py +++ b/quantflow_tests/test_ou.py @@ -1,6 +1,9 @@ +import numpy as np import pytest from quantflow.sp.ou import GammaOU, Vasicek +from quantflow.sp.poisson import CompoundPoissonProcess +from quantflow.utils.distributions import Exponential from quantflow_tests.utils import analytical_tests, characteristic_tests @@ -25,6 +28,58 @@ def test_sample(gamma_ou: GammaOU) -> None: assert paths.dt == 0.01 +def test_gamma_ou_sample_mean_unbiased() -> None: + # Use rate != stationary so the time-dependent mean is observable. + # Pre-fix the off-by-one in _advance under-biased the empirical mean by + # roughly kappa*dt; this test catches a ~2% bias at 20000 paths. + np.random.seed(0) + process = GammaOU( + rate=0.5, + kappa=2.0, + bdlp=CompoundPoissonProcess[Exponential]( + intensity=2.0, jumps=Exponential(decay=10.0) + ), + ) + paths = process.sample(20000, time_horizon=5.0, time_steps=500) + for t_idx, t in [(100, 1.0), (300, 3.0), (500, 5.0)]: + assert paths.data[t_idx].mean() == pytest.approx( + float(process.analytical_mean(t)), abs=0.005 + ) + + +def test_gamma_ou_transient_moments() -> None: + # rate != stationary mean so the time dependence is observable + process = GammaOU( + rate=0.5, + kappa=2.0, + bdlp=CompoundPoissonProcess[Exponential]( + intensity=2.0, jumps=Exponential(decay=10.0) + ), + ) + analytical_tests(process) + # also check the limit converges to the stationary moments + assert process.analytical_mean(100.0) == pytest.approx(0.2, abs=1e-6) + assert process.analytical_variance(100.0) == pytest.approx(0.02, abs=1e-6) + + +def test_gamma_ou_integrated_log_laplace() -> None: + process = GammaOU.create(rate=1.0, decay=10.0, kappa=1.0) + t = 1.0 + assert process.integrated_log_laplace(t, 0.0) == pytest.approx(0.0) + u = np.array([0.5, 1.0, 2.0]) + kappa = process.kappa + beta = process.beta + lam = process.intensity + f = (1 - np.exp(-kappa * t)) / kappa + expected = -u * process.rate * f + (lam / (u + beta * kappa)) * ( + beta * kappa * np.log((beta * kappa + u * f * kappa) / (beta * kappa)) + - u * kappa * t + ) + np.testing.assert_allclose( + process.integrated_log_laplace(t, u), expected, rtol=1e-12 + ) + + def test_vasicek(vasicek: Vasicek) -> None: m = vasicek.marginal(10) characteristic_tests(m) @@ -33,3 +88,12 @@ def test_vasicek(vasicek: Vasicek) -> None: assert m.variance() == pytest.approx(0.1) assert m.mean_from_characteristic() == pytest.approx(1.0, 1e-3) assert m.std_from_characteristic() == pytest.approx(m.std(), 1e-3) + + +def test_vasicek_negative() -> None: + # Vasicek admits negative initial values and mean levels + process = Vasicek(rate=-0.5, kappa=2.0, theta=-0.2) + m = process.marginal(1.0) + expected_mean = -0.5 * np.exp(-2.0) + (-0.2) * (1 - np.exp(-2.0)) + assert m.mean() == pytest.approx(expected_mean) + assert process.domain_range().lb == -np.inf diff --git a/quantflow_tests/test_poisson.py b/quantflow_tests/test_poisson.py index 9b637d6e..6dff6814 100644 --- a/quantflow_tests/test_poisson.py +++ b/quantflow_tests/test_poisson.py @@ -87,6 +87,13 @@ def test_dsp_pdf(dsp: DSP): np.testing.assert_almost_equal(pdf2[:32], pdf1) +def test_dsp_mean_scales_with_horizon(dsp: DSP): + # default CIR is stationary at rate=theta=1, so E[N_t] = t + for t in (0.5, 1.0, 2.0): + m = dsp.marginal(t) + assert m.mean_from_characteristic() == pytest.approx(t, rel=1e-4) + + def test_compound_create_double_exponential(): poi = CompoundPoissonProcess.create(DoubleExponential, jump_intensity=20, vol=0.5) assert poi.intensity == 20