From 0341bac01c2f32bf8f4a0ea6c40ee1c292fff2ea Mon Sep 17 00:00:00 2001 From: Luca Date: Mon, 4 May 2026 14:34:12 +0100 Subject: [PATCH] Fix Wiener typo --- app/hurst.py | 46 +++++++++---------- docs/api/index.md | 2 +- docs/api/sp/index.md | 2 +- docs/api/sp/weiner.md | 3 -- docs/api/sp/wiener.md | 3 ++ docs/examples/fft.py | 20 ++++---- ..._pricer.py => wiener_volatility_pricer.py} | 6 +-- docs/theory/convexity_correction.md | 4 +- docs/theory/inversion.md | 10 ++-- docs/theory/levy.md | 2 +- docs/tutorials/option_pricing.md | 4 +- mkdocs.yml | 2 +- notebooks/models/{weiner.md => wiener.md} | 12 ++--- quantflow/sp/base.py | 5 +- quantflow/sp/cir.py | 24 ++++++---- quantflow/sp/copula.py | 16 +++---- quantflow/sp/jump_diffusion.py | 10 ++-- quantflow/sp/ou.py | 6 +-- quantflow/sp/poisson.py | 2 +- quantflow/sp/{weiner.py => wiener.py} | 2 +- quantflow/utils/marginal.py | 7 ++- quantflow_tests/test_cir.py | 16 ++++++- quantflow_tests/test_copula.py | 36 +++++++++++++-- quantflow_tests/test_jump_diffusion.py | 29 ++++++++++++ quantflow_tests/test_ohlc.py | 4 +- quantflow_tests/test_options_pricer.py | 6 +-- quantflow_tests/test_ou.py | 11 +++++ quantflow_tests/test_poisson.py | 10 ++++ quantflow_tests/test_strategies.py | 4 +- .../{test_weiner.py => test_wiener.py} | 26 +++++------ 30 files changed, 213 insertions(+), 117 deletions(-) delete mode 100644 docs/api/sp/weiner.md create mode 100644 docs/api/sp/wiener.md rename docs/examples/{weiner_volatility_pricer.py => wiener_volatility_pricer.py} (71%) rename notebooks/models/{weiner.md => wiener.md} (81%) rename quantflow/sp/{weiner.py => wiener.py} (97%) rename quantflow_tests/{test_weiner.py => test_wiener.py} (68%) diff --git a/app/hurst.py b/app/hurst.py index 93c88ffc..3b296b87 100644 --- a/app/hurst.py +++ b/app/hurst.py @@ -40,9 +40,9 @@ def _(mo): We want to construct a mechanism to estimate the Hurst exponent via OHLC data because it is widely available from data providers and easily constructed as an online signal during trading. - In order to evaluate results against known solutions, we consider the Weiner process as generator of timeseries. + In order to evaluate results against known solutions, we consider the Wiener process as generator of timeseries. - We use the **WeinerProcess** from the stochastic process library and sample one path over a time horizon of 1 (day) with a time step every second. + We use the **WienerProcess** from the stochastic process library and sample one path over a time horizon of 1 (day) with a time step every second. """) return @@ -56,24 +56,24 @@ def _(mo): @app.cell def _(regenerate_btn): - from quantflow.sp.weiner import WeinerProcess + from quantflow.sp.wiener import WienerProcess from quantflow.utils import plot from quantflow.utils.dates import start_of_day regenerate_btn - weiner = WeinerProcess(sigma=2.0) - weiner_paths = weiner.sample(n=1, time_horizon=1, time_steps=24*60*60) - weiner_df = weiner_paths.as_datetime_df(start=start_of_day(), unit="d").reset_index() + wiener = WienerProcess(sigma=2.0) + wiener_paths = wiener.sample(n=1, time_horizon=1, time_steps=24*60*60) + wiener_df = wiener_paths.as_datetime_df(start=start_of_day(), unit="d").reset_index() plot.plot_lines( - weiner_df, - x=weiner_df.columns[0], - y=weiner_df.columns[1], - title="Weiner Process Path", - labels={"value": "Value", "variable": "Path", weiner_df.columns[0]: "Date"}, + wiener_df, + x=wiener_df.columns[0], + y=wiener_df.columns[1], + title="Wiener Process Path", + labels={"value": "Value", "variable": "Path", wiener_df.columns[0]: "Date"}, ) - return plot, start_of_day, weiner_df, weiner_paths + return plot, start_of_day, wiener_df, wiener_paths @app.cell @@ -90,14 +90,14 @@ def _(mo): ### Realized Variance At this point we estimate the standard deviation using the **realized variance** along the path (we use the **scaled** flag so that the standard deviation is scaled by the square-root of time step, in this way it removes the dependency on the time step size). - The value should be close to the **sigma** of the WeinerProcess defined above. + The value should be close to the **sigma** of the WienerProcess defined above. """) return @app.cell -def _(weiner_paths): - float(weiner_paths.paths_std(scaled=True)[0]) +def _(wiener_paths): + float(wiener_paths.paths_std(scaled=True)[0]) return @@ -110,15 +110,15 @@ def _(mo): @app.cell -def _(weiner_paths): - weiner_paths.hurst_exponent() +def _(wiener_paths): + wiener_paths.hurst_exponent() return @app.cell def _(mo): mo.md(r""" - As expected, the Hurst exponent should be close to 0.5, since we have calculated the exponent from the paths of a Weiner process. + As expected, the Hurst exponent should be close to 0.5, since we have calculated the exponent from the paths of a Wiener process. """) return @@ -143,7 +143,7 @@ def _(mo): @app.cell -def _(weiner_df): +def _(wiener_df): import pandas as pd import polars as pl import math @@ -167,7 +167,7 @@ def rstd(pdf: pl.Series, range_seconds: float) -> float: results = [] for period in ("10s", "20s", "30s", "1m", "2m", "3m", "5m", "10m", "30m"): ohlc = template.model_copy(update=dict(period=period)) - rf = ohlc(weiner_df) + rf = ohlc(wiener_df) ts = pd.to_timedelta(period).to_pytimedelta().total_seconds() data = dict(period=period) for name in ("pk", "gk", "rs"): @@ -255,15 +255,15 @@ def ohlc_hurst_exponent( @app.cell -def _(ohlc_hurst_exponent, weiner_df): - ohlc_hurst_exponent(weiner_df, series=["0"]) +def _(ohlc_hurst_exponent, wiener_df): + ohlc_hurst_exponent(wiener_df, series=["0"]) return @app.cell def _(mo): mo.md(r""" - The Hurst exponent should be close to 0.5, since we have calculated the exponent from the paths of a Weiner process. But it is not exactly 0.5 because the range-based estimators are not the same as the realized variance. Interestingly, the Parkinson estimator gives a Hurst exponent closer to 0.5 than the Garman-Klass and Rogers-Satchell estimators. + The Hurst exponent should be close to 0.5, since we have calculated the exponent from the paths of a Wiener process. But it is not exactly 0.5 because the range-based estimators are not the same as the realized variance. Interestingly, the Parkinson estimator gives a Hurst exponent closer to 0.5 than the Garman-Klass and Rogers-Satchell estimators. ## Mean Reverting Time Series diff --git a/docs/api/index.md b/docs/api/index.md index 00861a66..152aa9e6 100644 --- a/docs/api/index.md +++ b/docs/api/index.md @@ -37,7 +37,7 @@ Continuous-time stochastic processes used as underlying models for option pricin | Module | Description | |---|---| -| [Weiner Process](sp/weiner.md) | Geometric Brownian motion (constant volatility) | +| [Wiener Process](sp/wiener.md) | Geometric Brownian motion (constant volatility) | | [Heston Model](sp/heston.md) | Stochastic volatility with optional jump component (HestonJ) | | [Jump Diffusion](sp/jump_diffusion.md) | Compound Poisson jump processes | | [CIR Process](sp/cir.md) | Cox-Ingersoll-Ross mean-reverting process | diff --git a/docs/api/sp/index.md b/docs/api/sp/index.md index 8af83414..1d72867f 100644 --- a/docs/api/sp/index.md +++ b/docs/api/sp/index.md @@ -8,7 +8,7 @@ This page gives an overview of all stochastic processes available in the library | Process | Description | |---|---| -| [WeinerProcess][quantflow.sp.weiner.WeinerProcess] | Standard Brownian motion | +| [WienerProcess][quantflow.sp.wiener.WienerProcess] | Standard Brownian motion | ### Mean-reverting (intensity) diff --git a/docs/api/sp/weiner.md b/docs/api/sp/weiner.md deleted file mode 100644 index ff63446f..00000000 --- a/docs/api/sp/weiner.md +++ /dev/null @@ -1,3 +0,0 @@ -# Weiner process - -::: quantflow.sp.weiner.WeinerProcess diff --git a/docs/api/sp/wiener.md b/docs/api/sp/wiener.md new file mode 100644 index 00000000..43508857 --- /dev/null +++ b/docs/api/sp/wiener.md @@ -0,0 +1,3 @@ +# Wiener process + +::: quantflow.sp.wiener.WienerProcess diff --git a/docs/examples/fft.py b/docs/examples/fft.py index b0587fb4..c0c0743c 100644 --- a/docs/examples/fft.py +++ b/docs/examples/fft.py @@ -1,22 +1,22 @@ from docs.examples._utils import assets_path -from quantflow.sp.weiner import WeinerProcess +from quantflow.sp.wiener import WienerProcess from quantflow.utils import plot -p = WeinerProcess(sigma=0.5) +p = WienerProcess(sigma=0.5) m = p.marginal(0.2) fig = plot.plot_characteristic(m) -fig.update_layout(title="Weiner Process Characteristic Function") -fig.write_image(assets_path("weiner_characteristic.png")) +fig.update_layout(title="Wiener Process Characteristic Function") +fig.write_image(assets_path("wiener_characteristic.png")) fig = plot.plot_marginal_pdf(m, n=128, use_fft=True, max_frequency=20) -fig.update_layout(title="Weiner Process PDF via FFT with n=128") -fig.write_image(assets_path("weiner_fft_128.png")) +fig.update_layout(title="Wiener Process PDF via FFT with n=128") +fig.write_image(assets_path("wiener_fft_128.png")) fig = plot.plot_marginal_pdf(m, n=128 * 8, use_fft=True, max_frequency=8 * 20) -fig.update_layout(title="Weiner Process PDF via FFT with n=1024") -fig.write_image(assets_path("weiner_fft_1024.png")) +fig.update_layout(title="Wiener Process PDF via FFT with n=1024") +fig.write_image(assets_path("wiener_fft_1024.png")) fig = plot.plot_marginal_pdf(m, 64) -fig.update_layout(title="Weiner Process PDF via FRFT with n=64") -fig.write_image(assets_path("weiner_64.png")) +fig.update_layout(title="Wiener Process PDF via FRFT with n=64") +fig.write_image(assets_path("wiener_64.png")) diff --git a/docs/examples/weiner_volatility_pricer.py b/docs/examples/wiener_volatility_pricer.py similarity index 71% rename from docs/examples/weiner_volatility_pricer.py rename to docs/examples/wiener_volatility_pricer.py index eae572f5..472ad0ca 100644 --- a/docs/examples/weiner_volatility_pricer.py +++ b/docs/examples/wiener_volatility_pricer.py @@ -1,10 +1,10 @@ from quantflow.options.inputs import OptionType from quantflow.options.pricer import OptionPricer -from quantflow.sp.weiner import WeinerProcess +from quantflow.sp.wiener import WienerProcess -# Weiner process with constant volatility +# Wiener process with constant volatility # This produces the same sensitivities as the Black-Scholes model -pricer = OptionPricer(model=WeinerProcess(sigma=0.3)) +pricer = OptionPricer(model=WienerProcess(sigma=0.3)) # Price an ATM call option at time to maturity 1.0 price = pricer.price( diff --git a/docs/theory/convexity_correction.md b/docs/theory/convexity_correction.md index 4eb4a7ea..80454276 100644 --- a/docs/theory/convexity_correction.md +++ b/docs/theory/convexity_correction.md @@ -94,8 +94,8 @@ with time, reflecting the fact that geometric Brownian motion drifts upward on a (the asset price is log-normally distributed with positive variance). ```python -from quantflow.sp.weiner import WeinerProcess -pr = WeinerProcess(sigma=0.5) +from quantflow.sp.wiener import WienerProcess +pr = WienerProcess(sigma=0.5) -pr.characteristic_exponent(1, complex(0,-1)) # c_t at t=1 ``` diff --git a/docs/theory/inversion.md b/docs/theory/inversion.md index 8849aab3..acceceb1 100644 --- a/docs/theory/inversion.md +++ b/docs/theory/inversion.md @@ -56,17 +56,17 @@ which means $\delta_u$ and $\delta_x$ cannot be chosen independently — they ar \delta_x = \frac{2\pi}{N \delta_u} \end{equation} -As an example, let us invert the characteristic function of the Weiner process, which yields the standard normal distribution. +As an example, let us invert the characteristic function of the Wiener process, which yields the standard normal distribution. ```python --8<-- "docs/examples/fft.py" ``` -![Weiner Characteristic Function](../assets/examples/weiner_characteristic.png) +![Wiener Characteristic Function](../assets/examples/wiener_characteristic.png) -![Weiner FFT 128](../assets/examples/weiner_fft_128.png) +![Wiener FFT 128](../assets/examples/wiener_fft_128.png) -![Weiner FFT 1024](../assets/examples/weiner_fft_1024.png) +![Wiener FFT 1024](../assets/examples/wiener_fft_1024.png) **Note** the amount of unnecessary discretization points in the frequency domain (the characteristic function is zero after 15 or so). However the space domain is poorly represented because of the FFT constraints (we have a relatively small number of points where it matters, around zero). @@ -83,7 +83,7 @@ z &= \left(\left[e^{i j^2 \zeta/2}\right]_{j=0}^{N-1}, \left[e^{i\left(N-j\right We can now reduce the number of points needed for the discretization and achieve higher accuracy by properly selecting the domain discretization independently. -![Weiner FRFT 64](../assets/examples/weiner_64.png) +![Wiener FRFT 64](../assets/examples/wiener_64.png) Since one N-point FRFT will invoke three 2N-point FFT procedures, the number of operations will be approximately $6N\log{N}$ compared to $N\log{N}$ for the FFT. However, we can use fewer points as demonstrated and be more robust in delivering results. diff --git a/docs/theory/levy.md b/docs/theory/levy.md index 702c0663..bd9d8ce5 100644 --- a/docs/theory/levy.md +++ b/docs/theory/levy.md @@ -23,7 +23,7 @@ The independence and stationarity of the increments of the Lévy process imply t where the characteristic exponent $\phi_{x_1,u}$ is given by the [Lévy-Khintchine formula](https://en.wikipedia.org/wiki/L%C3%A9vy_process). There are several Lévy processes in the literature, including the [Poisson process](../api/sp/poisson.md), -the compound Poisson process, and the [Brownian motion](../api/sp/weiner.md). +the compound Poisson process, and the [Brownian motion](../api/sp/wiener.md). ## Time-Changed Lévy Processes diff --git a/docs/tutorials/option_pricing.md b/docs/tutorials/option_pricing.md index 30e224f3..fc0dab3d 100644 --- a/docs/tutorials/option_pricing.md +++ b/docs/tutorials/option_pricing.md @@ -13,11 +13,11 @@ Model sensitivities such as delta and gamma are calculated using the model dynam The first example shows how to price an option using the Black-Scholes model and validate the results against the analytical Black formula. The implied volatility should be the same as the model volatility, and the deltas and gammas should be the same as well (within numerical precision). ```python ---8<-- "docs/examples/weiner_volatility_pricer.py" +--8<-- "docs/examples/wiener_volatility_pricer.py" ``` ```json ---8<-- "docs/examples/output/weiner_volatility_pricer.out" +--8<-- "docs/examples/output/wiener_volatility_pricer.out" ``` diff --git a/mkdocs.yml b/mkdocs.yml index 2260d10b..526614a2 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -85,7 +85,7 @@ nav: - Jump Diffusion: api/sp/jump_diffusion.md - Ornstein-Uhlenbeck: api/sp/ou.md - Poisson Process: api/sp/poisson.md - - Weiner Process: api/sp/weiner.md + - Wiener Process: api/sp/wiener.md - Copulas: api/sp/copula.md - Technical Analysis: - api/ta/index.md diff --git a/notebooks/models/weiner.md b/notebooks/models/wiener.md similarity index 81% rename from notebooks/models/weiner.md rename to notebooks/models/wiener.md index 073eb66e..4adf68e5 100644 --- a/notebooks/models/weiner.md +++ b/notebooks/models/wiener.md @@ -11,9 +11,9 @@ kernelspec: name: python3 --- -# Weiner Process +# Wiener Process -In this document, we use the term Weiner process $w_t$ to indicate a Brownian motion with standard deviation given by the parameter $\sigma$; that is to say, the one-dimensional Weiner process is defined as: +In this document, we use the term Wiener process $w_t$ to indicate a Brownian motion with standard deviation given by the parameter $\sigma$; that is to say, the one-dimensional Wiener process is defined as: 1. $w_t$ is Lévy process 2. $d w_t = w_{t+dt}-w_t \sim N\left(0, \sigma dt\right)$ where $N$ is the normal distribution @@ -24,9 +24,9 @@ The [](characteristic-exponent) of $w$ is \end{equation} ```{code-cell} -from quantflow.sp.weiner import WeinerProcess +from quantflow.sp.wiener import WienerProcess -pr = WeinerProcess(sigma=0.5) +pr = WienerProcess(sigma=0.5) pr ``` @@ -47,8 +47,8 @@ plot.plot_marginal_pdf(m, 128) ```{code-cell} from quantflow.options.pricer import OptionPricer -from quantflow.sp.weiner import WeinerProcess -pricer = OptionPricer(WeinerProcess(sigma=0.2)) +from quantflow.sp.wiener import WienerProcess +pricer = OptionPricer(WienerProcess(sigma=0.2)) pricer ``` diff --git a/quantflow/sp/base.py b/quantflow/sp/base.py index 03c6ef9c..c46a1e32 100755 --- a/quantflow/sp/base.py +++ b/quantflow/sp/base.py @@ -178,7 +178,10 @@ def call_option_carr_madan_alpha(self) -> float: The choice of alpha is crucial for the numerical stability of the Carr-Madan formula. A common choice is to set alpha to a value that ensures the integrand - decays sufficiently fast at high frequencies. + decays sufficiently fast at high frequencies. The maturity-dependent heuristic + below (high alpha for short maturities, decaying to a 0.5 floor for long ones) + has been found to work well in practice across the processes in this library; + override on a per-process basis if needed. """ return max(8 * np.max(np.exp(-2 * self.t)), 0.5) diff --git a/quantflow/sp/cir.py b/quantflow/sp/cir.py index 900a9245..fd51be04 100755 --- a/quantflow/sp/cir.py +++ b/quantflow/sp/cir.py @@ -11,10 +11,10 @@ from .base import IntensityProcess -class SamplingAlgorithm(str, enum.Enum): - euler = "euler" - milstein = "milstein" - implicit = "implicit" +class SamplingAlgorithm(enum.StrEnum): + EULER = enum.auto() + MILSTEIN = enum.auto() + IMPLICIT = enum.auto() class CIR(IntensityProcess): @@ -36,7 +36,7 @@ class CIR(IntensityProcess): sigma: float = Field(default=1.0, gt=0, description=r"Volatility $\sigma$") theta: float = Field(default=1.0, gt=0, description=r"Mean rate $\theta$") sample_algo: SamplingAlgorithm = Field( - default=SamplingAlgorithm.implicit, description="Sampling algorithm" + default=SamplingAlgorithm.IMPLICIT, description="Sampling algorithm" ) @property @@ -64,14 +64,14 @@ def sample( def sample_from_draws(self, paths: Paths, *args: Paths) -> Paths: match self.sample_algo: - case SamplingAlgorithm.euler: + case SamplingAlgorithm.EULER: return self.sample_euler(paths) - case SamplingAlgorithm.milstein: + case SamplingAlgorithm.MILSTEIN: return self.sample_euler(paths, 0.25) - case SamplingAlgorithm.implicit: + case SamplingAlgorithm.IMPLICIT: return self.sample_implicit(paths) - def sample_euler(self, draws: Paths, ic: float = 0.0) -> Paths: + def sample_euler(self, draws: Paths, milstein_coef: float = 0.0) -> Paths: kappa = self.kappa theta = self.theta dt = draws.dt @@ -83,7 +83,11 @@ def sample_euler(self, draws: Paths, ic: float = 0.0) -> Paths: w = sdt * draws.data[t, :] x = paths[t, :] xplus = np.clip(x, 0, None) - dx = kappa * (theta - xplus) * dt + np.sqrt(xplus) * w + ic * (w * w - sdt2) + dx = ( + kappa * (theta - xplus) * dt + + np.sqrt(xplus) * w + + milstein_coef * (w * w - sdt2) + ) paths[t + 1, :] = x + dx return Paths(t=draws.t, data=paths) diff --git a/quantflow/sp/copula.py b/quantflow/sp/copula.py index e654e587..c50f6fce 100644 --- a/quantflow/sp/copula.py +++ b/quantflow/sp/copula.py @@ -1,12 +1,10 @@ from abc import ABC, abstractmethod -from decimal import Decimal from math import isclose import numpy as np from pydantic import BaseModel, Field from quantflow.utils.functions import debye -from quantflow.utils.numbers import ZERO from quantflow.utils.types import FloatArray, FloatArrayLike @@ -75,10 +73,10 @@ class FrankCopula(Copula): \end{equation} """ - kappa: Decimal = Field(default=ZERO, description="Frank copula parameter") + kappa: float = Field(default=0.0, description="Frank copula parameter") def __call__(self, u: FloatArrayLike, v: FloatArrayLike) -> FloatArrayLike: - k = float(self.kappa) + k = self.kappa if isclose(k, 0.0): return u * v eu = np.exp(-k * u) @@ -88,20 +86,20 @@ def __call__(self, u: FloatArrayLike, v: FloatArrayLike) -> FloatArrayLike: def tau(self) -> float: """Kendall's tau""" - k = float(self.kappa) + k = self.kappa if isclose(k, 0.0): return 0 return 1 + 4 * (debye(1, k) - 1) / k def rho(self) -> float: """Spearman's rho""" - k = float(self.kappa) + k = self.kappa if isclose(k, 0.0): return 0 return 1 - 12 * (debye(2, -k) - debye(1, -k)) / k def jacobian(self, u: FloatArrayLike, v: FloatArrayLike) -> FloatArray: - k = float(self.kappa) + k = self.kappa if isclose(k, 0.0): return np.array([v, u, v * 0]) eu = np.exp(-k * u) @@ -111,8 +109,6 @@ def jacobian(self, u: FloatArrayLike, v: FloatArrayLike) -> FloatArray: c = -np.log(1 + x) / k xx = x / (1 + x) du = eu * (ev - 1) / (e - 1) / (1 + x) - # du = eu * xx / (eu - 1) - dv = eu * (eu - 1) / (e - 1) / (1 + x) - # dv = ev * xx / (ev - 1) + dv = ev * (eu - 1) / (e - 1) / (1 + x) dk = (u * du + v * dv - e * xx / (e - 1) - c) / k return np.array([du, dv, dk]) diff --git a/quantflow/sp/jump_diffusion.py b/quantflow/sp/jump_diffusion.py index 9f480aa3..520f43e5 100644 --- a/quantflow/sp/jump_diffusion.py +++ b/quantflow/sp/jump_diffusion.py @@ -10,7 +10,7 @@ from ..utils.types import FloatArrayLike, Vector from .base import StochasticProcess1D from .poisson import CompoundPoissonProcess, D -from .weiner import WeinerProcess +from .wiener import WienerProcess class JumpDiffusion(StochasticProcess1D, Generic[D]): @@ -20,14 +20,14 @@ class JumpDiffusion(StochasticProcess1D, Generic[D]): dx_t = \sigma d w_t + d N_t \end{equation} - where $w_t$ is a [WeinerProcess][quantflow.sp.weiner.WeinerProcess] process + where $w_t$ is a [WienerProcess][quantflow.sp.wiener.WienerProcess] process with standard deviation $\sigma$ and $N_t$ is a [CompoundPoissonProcess][quantflow.sp.poisson.CompoundPoissonProcess] with intensity $\lambda$ and generic jump distribution $D$ """ - diffusion: WeinerProcess = Field( - default_factory=WeinerProcess, description="diffusion process" + diffusion: WienerProcess = Field( + default_factory=WienerProcess, description="diffusion process" ) jumps: CompoundPoissonProcess[D] = Field(description="The jump process") @@ -97,6 +97,6 @@ def create( jump_distribution_variance, jump_asymmetry ) return cls( - diffusion=WeinerProcess(sigma=np.sqrt(variance * (1 - jump_fraction))), + diffusion=WienerProcess(sigma=np.sqrt(variance * (1 - jump_fraction))), jumps=CompoundPoissonProcess(intensity=jump_intensity, jumps=jumps), ) diff --git a/quantflow/sp/ou.py b/quantflow/sp/ou.py index dde0fe97..cff1ca23 100755 --- a/quantflow/sp/ou.py +++ b/quantflow/sp/ou.py @@ -12,7 +12,7 @@ from ..utils.types import Float, FloatArrayLike, Vector from .base import IntensityProcess, StochasticProcess1D from .poisson import CompoundPoissonProcess, D -from .weiner import WeinerProcess +from .wiener import WienerProcess class Vasicek(StochasticProcess1D): @@ -43,8 +43,8 @@ class Vasicek(StochasticProcess1D): 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, + bdlp: WienerProcess = Field( + default_factory=WienerProcess, description="Background driving Wiener process", ) diff --git a/quantflow/sp/poisson.py b/quantflow/sp/poisson.py index 62b2be08..5ba20bc4 100755 --- a/quantflow/sp/poisson.py +++ b/quantflow/sp/poisson.py @@ -270,7 +270,7 @@ def cdf_from_characteristic( c = self.characteristic(frequency) a = 1 / np.pi result = [] - x = np.arange(n or 10) + x = np.arange(n) for m in x: d = np.sin(0.5 * frequency) d[0] = 1.0 diff --git a/quantflow/sp/weiner.py b/quantflow/sp/wiener.py similarity index 97% rename from quantflow/sp/weiner.py rename to quantflow/sp/wiener.py index 877a49d0..ed016c87 100644 --- a/quantflow/sp/weiner.py +++ b/quantflow/sp/wiener.py @@ -9,7 +9,7 @@ from .base import StochasticProcess1D -class WeinerProcess(StochasticProcess1D): +class WienerProcess(StochasticProcess1D): sigma: float = Field(default=1, ge=0, description="volatility") @property diff --git a/quantflow/utils/marginal.py b/quantflow/utils/marginal.py index d68c4e36..e72365b9 100644 --- a/quantflow/utils/marginal.py +++ b/quantflow/utils/marginal.py @@ -215,8 +215,11 @@ def call_option_carr_madan( def call_option_carr_madan_alpha(self) -> float: """Option alpha to use for Carr & Madan transform to ensure integrability - of the call option transform""" - return 2.0 + of the call option transform. + + Defaults to 1.5, the value suggested in the original Carr & Madan paper. + """ + return 1.5 def call_option_cos( self, diff --git a/quantflow_tests/test_cir.py b/quantflow_tests/test_cir.py index 5033771c..82f27437 100644 --- a/quantflow_tests/test_cir.py +++ b/quantflow_tests/test_cir.py @@ -6,12 +6,12 @@ @pytest.fixture def cir_neg() -> CIR: - return CIR(kappa=1, sigma=2, sample_algo=SamplingAlgorithm.euler) + return CIR(kappa=1, sigma=2, sample_algo=SamplingAlgorithm.EULER) @pytest.fixture def cir() -> CIR: - return CIR(kappa=1, sigma=1.2, sample_algo=SamplingAlgorithm.euler) + return CIR(kappa=1, sigma=1.2, sample_algo=SamplingAlgorithm.EULER) def test_feller_condition(cir: CIR, cir_neg: CIR) -> None: @@ -43,3 +43,15 @@ def test_cir_pdf(cir: CIR): m = cir.marginal(1) pdf = m.pdf_from_characteristic(128, max_frequency=20) np.testing.assert_array_almost_equal(pdf.y, m.pdf(pdf.x), 1e-1) + + +def test_cir_pdf_neg(cir_neg: CIR): + # q = 2 kappa theta / sigma^2 - 1 = -0.5 < 0 (Feller violated) + assert cir_neg.is_positive is False + q = 2 * cir_neg.kappa * cir_neg.theta / cir_neg.sigma2 - 1 + assert q < 0 + m = cir_neg.marginal(1) + pdf = m.pdf_from_characteristic(128, max_frequency=20) + # skip x=0 where the analytical pdf diverges for q<0 + mask = pdf.x > 0 + np.testing.assert_array_almost_equal(pdf.y[mask], m.pdf(pdf.x[mask]), 1e-1) diff --git a/quantflow_tests/test_copula.py b/quantflow_tests/test_copula.py index 93966f6b..f770048e 100644 --- a/quantflow_tests/test_copula.py +++ b/quantflow_tests/test_copula.py @@ -1,4 +1,3 @@ -from decimal import Decimal from math import isclose import numpy as np @@ -14,8 +13,8 @@ def test_independent_copula(): def test_frank_copula(): - c = FrankCopula(kappa=Decimal("0.3")) - assert c.kappa == Decimal("0.3") + c = FrankCopula(kappa=0.3) + assert c.kappa == 0.3 assert c.tau() > 0 assert c.rho() < 0 assert c.jacobian(0.3, 0.4).shape == (3,) @@ -27,4 +26,33 @@ def test_frank_copula(): c = FrankCopula() assert isclose(c(11.0, 3.0), 33.0) - assert isclose(c(11.0, 3.0), 33.0) + + +def test_frank_copula_call(): + c = FrankCopula(kappa=2.5) + # boundaries: C(0, v) = 0, C(u, 0) = 0, C(1, v) = v, C(u, 1) = u + assert isclose(c(0.0, 0.4), 0.0, abs_tol=1e-12) + assert isclose(c(0.3, 0.0), 0.0, abs_tol=1e-12) + assert isclose(c(1.0, 0.4), 0.4) + assert isclose(c(0.3, 1.0), 0.3) + # symmetry: C(u, v) = C(v, u) + assert isclose(c(0.3, 0.7), c(0.7, 0.3)) + # bounded by Frechet: max(u+v-1, 0) <= C(u, v) <= min(u, v) + assert max(0.3 + 0.4 - 1, 0) <= c(0.3, 0.4) <= min(0.3, 0.4) + + +def test_frank_copula_jacobian(): + c = FrankCopula(kappa=2.5) + u, v = 0.3, 0.4 + h = 1e-6 + du, dv, dk = c.jacobian(u, v) + # finite differences vs analytical partials + du_num = (c(u + h, v) - c(u - h, v)) / (2 * h) + dv_num = (c(u, v + h) - c(u, v - h)) / (2 * h) + assert isclose(float(du), du_num, abs_tol=1e-7) + assert isclose(float(dv), dv_num, abs_tol=1e-7) + # dk: numerical via perturbing kappa + c_up = FrankCopula(kappa=2.5 + 1e-6) + c_dn = FrankCopula(kappa=2.5 - 1e-6) + dk_num = (c_up(u, v) - c_dn(u, v)) / 2e-6 + assert isclose(float(dk), dk_num, abs_tol=1e-6) diff --git a/quantflow_tests/test_jump_diffusion.py b/quantflow_tests/test_jump_diffusion.py index c6fb00c8..168ec4c8 100644 --- a/quantflow_tests/test_jump_diffusion.py +++ b/quantflow_tests/test_jump_diffusion.py @@ -1,4 +1,6 @@ +import numpy as np import pytest +from scipy.stats import norm, poisson from quantflow.sp.jump_diffusion import JumpDiffusion from quantflow.utils.distributions import Normal @@ -25,3 +27,30 @@ def test_sampling(merton: JumpDiffusion[Normal]) -> None: assert mean[0] == 0 std = paths.std() assert std[0] == 0 + + +def merton_analytical_pdf( + merton: JumpDiffusion[Normal], t: float, x: np.ndarray, max_n: int = 60 +) -> np.ndarray: + # Merton model: x_t = sigma*W_t + sum_{i=1}^{N_t} J_i with J_i ~ N(mu, sigma_J^2), + # density is the Poisson-weighted mixture of Gaussians. + sigma = merton.diffusion.sigma + lam = merton.jumps.intensity + mu_j = merton.jumps.jumps.mu + sigma_j = merton.jumps.jumps.sigma + pdf = np.zeros_like(x, dtype=float) + for n in range(max_n): + var = sigma * sigma * t + n * sigma_j * sigma_j + pdf += poisson.pmf(n, lam * t) * norm.pdf(x, loc=n * mu_j, scale=np.sqrt(var)) + return pdf + + +def test_pdf_vs_characteristic(merton: JumpDiffusion[Normal]) -> None: + t = 1.0 + m = merton.marginal(t) + pdf = m.pdf_from_characteristic(256) + analytical = merton_analytical_pdf(merton, t, np.asarray(pdf.x)) + # restrict to the bulk where the analytical density is non-negligible: + # the FFT-based reconstruction has tail aliasing where pdf ~ 0. + mask = analytical > 1e-3 + np.testing.assert_array_almost_equal(pdf.y[mask], analytical[mask], decimal=2) diff --git a/quantflow_tests/test_ohlc.py b/quantflow_tests/test_ohlc.py index d3701ea2..dd9c7e6c 100644 --- a/quantflow_tests/test_ohlc.py +++ b/quantflow_tests/test_ohlc.py @@ -1,4 +1,4 @@ -from quantflow.sp.weiner import WeinerProcess +from quantflow.sp.wiener import WienerProcess from quantflow.ta.ohlc import OHLC @@ -18,7 +18,7 @@ def test_ohlc() -> None: assert ohlc.rogers_satchell_variance is True assert ohlc.percent_variance is False # create a dataframe - path = WeinerProcess(sigma=0.5).sample(1, 1, 1000) + path = WienerProcess(sigma=0.5).sample(1, 1, 1000) df = path.as_datetime_df().reset_index() result = ohlc(df) assert result.shape == (145, 9) diff --git a/quantflow_tests/test_options_pricer.py b/quantflow_tests/test_options_pricer.py index 2762e915..08e2de60 100644 --- a/quantflow_tests/test_options_pricer.py +++ b/quantflow_tests/test_options_pricer.py @@ -4,7 +4,7 @@ from quantflow.options.pricer import OptionPricer, OptionType from quantflow.sp.heston import HestonJ -from quantflow.sp.weiner import WeinerProcess +from quantflow.sp.wiener import WienerProcess from quantflow.utils.distributions import DoubleExponential from quantflow_tests.utils import has_plotly @@ -41,9 +41,9 @@ def test_price_call(pricer: OptionPricer): @pytest.mark.parametrize("strike,forward", [(90, 100), (100, 100), (110, 100)]) -def test_weiner_matches_black(strike: int, forward: int) -> None: +def test_wiener_matches_black(strike: int, forward: int) -> None: sigma = 0.3 - pricer = OptionPricer(model=WeinerProcess(sigma=sigma)) + pricer = OptionPricer(model=WienerProcess(sigma=sigma)) price = pricer.price( option_type=OptionType.call, strike=strike, forward=forward, ttm=1.0 ) diff --git a/quantflow_tests/test_ou.py b/quantflow_tests/test_ou.py index 62a822fe..9cd76be9 100644 --- a/quantflow_tests/test_ou.py +++ b/quantflow_tests/test_ou.py @@ -1,5 +1,6 @@ import numpy as np import pytest +from scipy.integrate import cumulative_trapezoid from quantflow.sp.ou import GammaOU, Vasicek from quantflow.sp.poisson import CompoundPoissonProcess @@ -97,3 +98,13 @@ def test_vasicek_negative() -> None: 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 + + +def test_vasicek_pdf_cdf_parity(vasicek: Vasicek) -> None: + m = vasicek.marginal(1.0) + pdf = m.pdf_from_characteristic(128) + np.testing.assert_array_almost_equal(pdf.y, m.pdf(pdf.x), decimal=2) + # analytical CDF is consistent with analytical PDF via numerical integration + x = np.linspace(m.mean() - 6 * m.std(), m.mean() + 6 * m.std(), 4096) + cdf_from_pdf = cumulative_trapezoid(m.pdf(x), x, initial=0) + m.cdf(x[0]) + np.testing.assert_array_almost_equal(m.cdf(x), cdf_from_pdf, decimal=4) diff --git a/quantflow_tests/test_poisson.py b/quantflow_tests/test_poisson.py index 6dff6814..dab84a6f 100644 --- a/quantflow_tests/test_poisson.py +++ b/quantflow_tests/test_poisson.py @@ -4,6 +4,7 @@ import pytest from quantflow.sp.dsp import DSP +from quantflow.sp.ou import GammaOU from quantflow.sp.poisson import CompoundPoissonProcess, PoissonProcess from quantflow.utils.distributions import DoubleExponential, Exponential, Normal from quantflow_tests.utils import analytical_tests, characteristic_tests @@ -94,6 +95,15 @@ def test_dsp_mean_scales_with_horizon(dsp: DSP): assert m.mean_from_characteristic() == pytest.approx(t, rel=1e-4) +def test_dsp_with_gamma_ou_intensity(): + # GammaOU stationary mean = intensity * E[jump] = (rate * decay) / decay = rate + intensity = GammaOU.create(rate=1.0, decay=2.0, kappa=1.5) + dsp = DSP(intensity=intensity) + for t in (0.5, 1.0, 2.0): + m = dsp.marginal(t) + assert m.mean_from_characteristic() == pytest.approx(t, rel=1e-3) + + def test_compound_create_double_exponential(): poi = CompoundPoissonProcess.create(DoubleExponential, jump_intensity=20, vol=0.5) assert poi.intensity == 20 diff --git a/quantflow_tests/test_strategies.py b/quantflow_tests/test_strategies.py index d7b07b97..80552162 100644 --- a/quantflow_tests/test_strategies.py +++ b/quantflow_tests/test_strategies.py @@ -10,7 +10,7 @@ Straddle, Strangle, ) -from quantflow.sp.weiner import WeinerProcess +from quantflow.sp.wiener import WienerProcess REF_DATE = datetime(2024, 1, 1, tzinfo=timezone.utc) MATURITY = datetime(2025, 1, 1, tzinfo=timezone.utc) @@ -20,7 +20,7 @@ @pytest.fixture def pricer() -> OptionPricer: - return OptionPricer(model=WeinerProcess(sigma=0.3)) + return OptionPricer(model=WienerProcess(sigma=0.3)) def test_straddle_from_strike(pricer: OptionPricer) -> None: diff --git a/quantflow_tests/test_weiner.py b/quantflow_tests/test_wiener.py similarity index 68% rename from quantflow_tests/test_weiner.py rename to quantflow_tests/test_wiener.py index e789c34c..af9dc561 100644 --- a/quantflow_tests/test_weiner.py +++ b/quantflow_tests/test_wiener.py @@ -1,19 +1,19 @@ import numpy as np import pytest -from quantflow.sp.weiner import WeinerProcess +from quantflow.sp.wiener import WienerProcess from quantflow_tests.utils import characteristic_tests @pytest.fixture -def weiner() -> WeinerProcess: - return WeinerProcess(sigma=0.5) +def wiener() -> WienerProcess: + return WienerProcess(sigma=0.5) -def test_characteristic(weiner: WeinerProcess) -> None: - assert weiner.characteristic(1, 0) == 1 - assert weiner.convexity_correction(2) == 0.25 - marginal = weiner.marginal(1) +def test_characteristic(wiener: WienerProcess) -> None: + assert wiener.characteristic(1, 0) == 1 + assert wiener.convexity_correction(2) == 0.25 + marginal = wiener.marginal(1) characteristic_tests(marginal) assert marginal.mean() == 0 assert marginal.mean_from_characteristic() == 0 @@ -24,22 +24,22 @@ def test_characteristic(weiner: WeinerProcess) -> None: assert len(df.columns) == 3 -def test_sampling(weiner: WeinerProcess) -> None: - paths = weiner.sample(1000, time_horizon=1, time_steps=1000) +def test_sampling(wiener: WienerProcess) -> None: + paths = wiener.sample(1000, time_horizon=1, time_steps=1000) mean = paths.mean() assert mean[0] == 0 std = paths.std() assert std[0] == 0 -def test_support(weiner: WeinerProcess) -> None: - m = weiner.marginal(0.01) +def test_support(wiener: WienerProcess) -> None: + m = wiener.marginal(0.01) pdf = m.pdf_from_characteristic(32) assert len(pdf.x) == 32 -def test_fft_v_frft(weiner: WeinerProcess) -> None: - m = weiner.marginal(1) +def test_fft_v_frft(wiener: WienerProcess) -> None: + m = wiener.marginal(1) pdf1 = m.pdf_from_characteristic(128, max_frequency=10) pdf2 = m.pdf_from_characteristic(128, use_fft=True, max_frequency=200) y = np.interp(pdf1.x[10:-10], pdf2.x, pdf2.y)