diff --git a/Makefile b/Makefile index d3a7e66c..8e4dbd74 100644 --- a/Makefile +++ b/Makefile @@ -11,6 +11,10 @@ docs: ## build documentation @uv run ./dev/build-examples @uv run mkdocs build +.PHONY: docs-bib +docs-bib: ## Regenerate docs bibliography + @uv run ./docs/bib2md.py + .PHONY: docs-examples docs-examples: ## Regenerate docs examples @uv run ./dev/build-examples diff --git a/docs/api/options/pricer.md b/docs/api/options/pricer.md index 796736d8..3b90ad1c 100644 --- a/docs/api/options/pricer.md +++ b/docs/api/options/pricer.md @@ -5,6 +5,8 @@ different stochastic volatility models. ::: quantflow.options.pricer.OptionPricerBase +::: quantflow.options.pricer.OptionPricingMethod + ::: quantflow.options.pricer.OptionPricer ::: quantflow.options.pricer.MaturityPricer diff --git a/docs/api/sp/bns.md b/docs/api/sp/bns.md new file mode 100644 index 00000000..ddbe84ff --- /dev/null +++ b/docs/api/sp/bns.md @@ -0,0 +1,3 @@ +# Barndorff-Nielson & Shephard process + +::: quantflow.sp.bns.BNS diff --git a/docs/api/sp/index.md b/docs/api/sp/index.md index d10439fd..8af83414 100644 --- a/docs/api/sp/index.md +++ b/docs/api/sp/index.md @@ -48,3 +48,5 @@ This page gives an overview of all stochastic processes available in the library ::: quantflow.sp.base.StochasticProcess1D ::: quantflow.sp.base.IntensityProcess + +::: quantflow.sp.base.StochasticProcess1DMarginal diff --git a/docs/api/sp/ou.md b/docs/api/sp/ou.md index 761e50e9..01e3bb34 100644 --- a/docs/api/sp/ou.md +++ b/docs/api/sp/ou.md @@ -10,4 +10,6 @@ These are the classes that implement gaussian and non-gaussian ## Non-Gaussian OU process +::: quantflow.sp.ou.NGOU + ::: quantflow.sp.ou.GammaOU diff --git a/docs/api/utils/index.md b/docs/api/utils/index.md index e9f445d3..921026f5 100644 --- a/docs/api/utils/index.md +++ b/docs/api/utils/index.md @@ -1,3 +1,15 @@ # Utilities -This section contains utility functions and classes that are used throughout the library. They are not meant to be used directly by the user, but they can be useful for advanced users who want to extend the library or understand its inner workings. +This section contains utility functions and classes used throughout the library. +They are not meant to be used directly by the user, but can be useful for advanced +users who want to extend the library or understand its inner workings. + +## Classes + +| Page | Description | +|---|---| +| [Marginal 1D](marginal1d.md) | Abstract base class for 1D marginal distributions with Fourier-based option pricing (Carr-Madan and Lewis formulas) | +| [Distributions](distributions.md) | Parametric 1D distributions (e.g. Exponential) | +| [Bins](bins.md) | Histogram and event-density utilities | +| [Numbers](numbers.md) | Decimal number helpers | +| [Types](types.md) | Shared type aliases (FloatArray, etc.) | diff --git a/docs/bib2md.py b/docs/bib2md.py new file mode 100644 index 00000000..eaecbd87 --- /dev/null +++ b/docs/bib2md.py @@ -0,0 +1,195 @@ +#!/usr/bin/env python3 +"""Convert docs/references.bib to docs/bibliography.md. + +Usage: + uv run python docs/bib2md.py [--bib PATH] [--out PATH] + +The BibTeX file is the source of truth; run this script whenever it changes. +""" +from __future__ import annotations + +import argparse +import re +import sys +from pathlib import Path + +# --------------------------------------------------------------------------- +# Parser +# --------------------------------------------------------------------------- + +_ENTRY_RE = re.compile( + r"@(\w+)\s*\{\s*([\w][\w:/-]*)\s*,\s*(.*?)\n\}", + re.DOTALL, +) +_FIELD_RE = re.compile( + r"""(\w+)\s*=\s*(?:\{((?:[^{}]|\{[^{}]*\})*)\}|"([^"]*)")""", + re.DOTALL, +) + + +def _strip_braces(text: str) -> str: + """Remove single-level LaTeX capitalisation braces, e.g. {Lévy} -> Lévy.""" + return re.sub(r"\{([^{}]*)\}", r"\1", text) + + +def _clean(value: str) -> str: + text = _strip_braces(" ".join(value.split())) + # Convert LaTeX en-dash (--) to a plain hyphen + return text.replace("--", "-") + + +def parse_bib(text: str) -> list[dict[str, str]]: + entries = [] + for m in _ENTRY_RE.finditer(text): + entry: dict[str, str] = { + "type": m.group(1).lower(), + "key": m.group(2), + } + for fm in _FIELD_RE.finditer(m.group(3)): + raw = fm.group(2) if fm.group(2) is not None else fm.group(3) + entry[fm.group(1).lower()] = _clean(raw) + entries.append(entry) + return entries + + +# --------------------------------------------------------------------------- +# Author formatting +# --------------------------------------------------------------------------- + +def _format_single_author(author: str) -> str: + """Convert 'Lastname, Firstname' to 'Firstname Lastname'; leave others as-is.""" + author = author.strip() + # Only reformat if a comma is clearly separating last from first name + # (avoid splitting on commas inside name suffixes or freeform strings) + parts = author.split(",", 1) + if len(parts) == 2: + last, first = parts[0].strip(), parts[1].strip() + # Sanity: last name should be a single word (no spaces) + if first and " " not in last: + return f"{first} {last}" + return author + + +def format_authors(author_field: str) -> str: + # Normalise separators: "and" and "&" both split authors + normalised = re.sub(r"\s+and\s+", " & ", author_field, flags=re.IGNORECASE) + authors = [a.strip() for a in normalised.split("&") if a.strip()] + return ", ".join(_format_single_author(a) for a in authors) + + +# --------------------------------------------------------------------------- +# Entry formatter +# --------------------------------------------------------------------------- + +def _title_link(title: str, url: str) -> str: + if not url: + return title + return f'[{title}]({url}){{target="_blank" rel="noopener"}}' + + +def _journal_detail(entry: dict[str, str]) -> str: + """Build 'Journal, Vol(Num):Pages' string.""" + journal = entry.get("journal", "") + volume = entry.get("volume", "") + number = entry.get("number", "") + pages = entry.get("pages", "") + if not journal: + return "" + detail = journal + if volume: + vol_str = volume + if number: + vol_str += f"({number})" + if pages: + vol_str += f":{pages}" + detail += f", {vol_str}" + elif pages: + detail += f", {pages}" + return detail + + +def format_entry(entry: dict[str, str]) -> str: + key = entry["key"] + title = entry.get("title", "Unknown Title") + year = entry.get("year", "") + url = entry.get("url", "") + author = entry.get("author", "") + + # Prefer doi as url when no url field present + if not url and "doi" in entry: + doi = entry["doi"].lstrip("doi:") + url = f"https://doi.org/{doi}" + + author_str = format_authors(author) if author else "" + title_md = _title_link(title, url) + year_str = f"({year})" if year else "" + + # Prefix: "Author(s). (Year)" + prefix_parts = [] + if author_str: + prefix_parts.append(author_str.rstrip(".") + ".") + if year_str: + prefix_parts.append(year_str) + prefix = " ".join(prefix_parts) + + # Suffix depends on entry type + entry_type = entry.get("type", "article") + if entry_type == "article": + suffix = _journal_detail(entry) + elif entry_type == "book": + suffix = entry.get("publisher", "") + elif entry_type in ("mastersthesis", "phdthesis"): + suffix = entry.get("school", "") + else: + suffix = entry.get("publisher", entry.get("school", "")) + + body = f"{prefix} {title_md}" if prefix else title_md + if suffix: + body = f"{body}, {suffix}" + + return f"#### {key}\n\n{body}\n" + + +# --------------------------------------------------------------------------- +# Main +# --------------------------------------------------------------------------- + +def convert(bib_path: Path, out_path: Path) -> None: + text = bib_path.read_text(encoding="utf-8") + entries = parse_bib(text) + if not entries: + print(f"No entries found in {bib_path}", file=sys.stderr) + sys.exit(1) + + entries.sort(key=lambda e: e["key"].lower()) + + lines = ["# Bibliography\n", "\n---\n"] + for entry in entries: + lines.append("\n") + lines.append(format_entry(entry)) + + out_path.write_text("".join(lines), encoding="utf-8") + print(f"Wrote {len(entries)} entries to {out_path}") + + +def main() -> None: + docs = Path(__file__).parent + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument( + "--bib", + type=Path, + default=docs / "references.bib", + help="Input BibTeX file (default: docs/references.bib)", + ) + parser.add_argument( + "--out", + type=Path, + default=docs / "bibliography.md", + help="Output Markdown file (default: docs/bibliography.md)", + ) + args = parser.parse_args() + convert(args.bib, args.out) + + +if __name__ == "__main__": + main() diff --git a/docs/bibliography.md b/docs/bibliography.md index 42371196..32cfb136 100644 --- a/docs/bibliography.md +++ b/docs/bibliography.md @@ -1,23 +1,73 @@ # Bibliography +The raw BibTeX source for all entries is available in [references.bib](/references.bib). + --- +#### bertoin + +J. Bertoin. (1996) Lévy Processes, Cambridge University Press + +#### bns + +O.E. Barndorff-Nielsen, N. Shephard. (2001) [Non-Gaussian OU based models and some of their uses in financial economics](https://drive.google.com/file/d/11TnGl2ER_-QzEu_h7ZtGGass2rF8_u1Y/view){target="_blank" rel="noopener"}, Journal of the Royal Statistical Society, 63(2) + #### carr_madan -Peter Carr, Dilip Madan, [Option Valuation Using the Fast Fourier Transform](https://doi.org/10.21314/JCF.1999.043), Journal of Computational Finance, 2(4):61-73, 1999 +P. Carr, D. Madan. (1999) [Option valuation using the fast Fourier transform](http://faculty.baruch.cuny.edu/lwu/890/CarrMadan99.pdf){target="_blank" rel="noopener"}, Journal of Computational Finance, 3:463-520 #### carr_wu -Peter Carr, Liuren Wu, [Time-Changed Lévy Processes and Option Pricing](https://doi.org/10.1016/S0304-405X(03)00171-5), Journal of Financial Economics, 71(1):113-141, 2004 +P. Carr, L. Wu. (2002) [Time-changed Lévy processes and option pricing](https://engineering.nyu.edu/sites/default/files/2019-03/Carr-time-changed-levy-processes-option-pricing.pdf){target="_blank" rel="noopener"}, Journal of Financial Economics, 7:113-141 + +#### cgmy + +Carr P., Geman H., Madan D.B., Yor M. (2003) [Stochastic Volatility for Lévy processes](https://engineering.nyu.edu/sites/default/files/2019-03/Carr-stochastic-volatility-levy-processes.pdf){target="_blank" rel="noopener"}, Mathematical Finance, 13(3) #### chourdakis -Kyriakos Chourdakis, [Option Pricing Using the Fractional FFT](https://doi.org/10.21314/JCF.2005.137),Journal of Computational Finance, 8(2):1-18, 2005 +K. Chourdakis. (2004) [Option Pricing Using the Fractional FFT](https://citeseerx.ist.psu.edu/document?repid=rep1&type=pdf&doi=6bdf4696312d37427eda2740137650c09deacda7){target="_blank" rel="noopener"}, Journal of Computational Finance, 8:1-18 + +#### cos + +F. Fang, K. Oosterlee. (2008) [A Novel Pricing Method for European Options Based on Fourier-Cosine Series Expansions](https://mpra.ub.uni-muenchen.de/8914/4/MPRA_paper_8914.pdf){target="_blank" rel="noopener"} + +#### dsp + +Unknown. (2017) [Doubly Stochastic Poisson Processes with Affine Intensities](https://editorialexpress.com/cgi-bin/conference/download.cgi?db_name=sbe35&paper_id=179){target="_blank" rel="noopener"}, internet + +#### ekf + +Wang X., He X., Zhao Y., Zuo Z. (2017) [Parameter Estimations of Heston Model Based on Consistent Extended Kalman Filter](https://www.sciencedirect.com/science/article/pii/S2405896317324758){target="_blank" rel="noopener"}, internet + +#### gamma-ou + +P. Sabino, C. Petroni. (2021) [Gamma Related Ornstein-Uhlenbeck Processes and their Simulation](https://doi.org/10.1080/00949655.2020.1842408){target="_blank" rel="noopener"}, Journal of Statistical Computation and Simulation, 91(6) + +#### heston-calibration + +Milan Mrázek, Jan Pospíšil. (2017) [Calibration and simulation of Heston model](https://doi.org/10.1515/math-2017-0058){target="_blank" rel="noopener"}, Open Mathematics, 15(1):679-704 + +#### heston-simulation + +Leif B.G. Andersen. (2008) [Efficient Simulation of the Heston Stochastic Volatility Model](https://papers.ssrn.com/sol3/papers.cfm?abstract_id=946405){target="_blank" rel="noopener"}, Journal of Computational Finance, 11(3) + +#### lee_option + +R. W. Lee. (2004) [Option Pricing by Transform Methods: Extensions, Unification, and Error Control](https://www.math.uchicago.edu/~rl/dft.pdf){target="_blank" rel="noopener"} + +#### lewis + +A. L. Lewis. (2001) [A Simple Option Formula for General Jump-Diffusion and other Exponential Lévy Processes](https://drive.google.com/file/d/1JiYfyOU7lUKHrMqTRPbBCdasctwYDc4Z/view){target="_blank" rel="noopener"} #### molnar -Peter Molnar, [Volatility modeling and forecasting: utilization of realized volatility, implied volatility and the highest and lowest price of the day](https://drive.google.com/file/d/1zCU1OZyrKQLpxaypPv9U5UPbReBDXcMf/view), Master's thesis, University of Economics in Prague, 2020 +Peter Molnar. (2020) [Volatility modeling and forecasting: utilization of realized volatility, implied volatility and the highest and lowest price of the day](https://drive.google.com/file/d/1zCU1OZyrKQLpxaypPv9U5UPbReBDXcMf/view){target="_blank" rel="noopener"}, University of Economics in Prague + +#### saez + +G. K. G. Saez. (2014) [Fourier Transform Methods for Option Pricing: An Application to extended Heston-type Models](https://www.uv.es/bfc/TFM2014/008-014.pdf){target="_blank" rel="noopener"}, Universidad del Pais Vasco -#### Gauthier +#### ukf -Gauthier Godin & Legros, [Deep Implied Volatility Factor Models for Stock Options](https://drive.google.com/file/d/1Rjypn1IqnhpiZz08s0hxl5ISQDC8KeWk/view?usp=sharing), 2025 +Merwe. (2014) [The Unscented Kalman Filter for Nonlinear Estimation](https://groups.seas.harvard.edu/courses/cs281/papers/unscented.pdf){target="_blank" rel="noopener"}, internet diff --git a/docs/examples/pricing_method_comparison.py b/docs/examples/pricing_method_comparison.py new file mode 100644 index 00000000..85d8f736 --- /dev/null +++ b/docs/examples/pricing_method_comparison.py @@ -0,0 +1,248 @@ +"""Compare Carr-Madan, Lewis and COS option pricing methods for accuracy and speed.""" + +import time + +import numpy as np +import plotly.graph_objects as go +from plotly.subplots import make_subplots +from pydantic import BaseModel, Field + +from docs.examples._utils import assets_path +from quantflow.options.bs import implied_black_volatility +from quantflow.sp.base import StochasticProcess1D +from quantflow.sp.heston import Heston +from quantflow.utils.marginal import OptionPricingMethod, OptionPricingResult + + +class ChartProps(BaseModel): + color: str = Field(default="#1f77b4", description="Line color for the chart") + dash: str = Field( + default="solid", + description="Line dash style for the chart (e.g., 'solid', 'dash', 'dot')", + ) + + +class PricingMethodComparison(BaseModel): + model: StochasticProcess1D = Field( + description="Stochastic process model to compare" + ) + ttms: tuple[float, ...] = Field( + default=(1.0, 0.25, 0.02), + description="Time to maturities to compare", + ) + ns: tuple[int, ...] = Field( + default=(32, 64, 128, 256, 512), + description="Discretization points to compare", + ) + timing_reps: int = Field( + default=20, + description="Number of repetitions for timing measurements", + ) + max_moneyness: float = Field( + default=1.5, + description="Maximum time-adjusted moneyness for option pricing", + ) + ref_n: int = Field( + default=8192, + description=( + "Number of discretization points to use for the reference Lewis price" + ), + ) + max_iv_error: float = Field( + default=0.1, + description="Implied vol errors above this value are clipped in the error plot", + ) + charts: dict[OptionPricingMethod, ChartProps] = Field( + default_factory=lambda: { + OptionPricingMethod.CARR_MADAN: ChartProps(color="#1f77b4", dash="solid"), + OptionPricingMethod.LEWIS: ChartProps(color="#ff7f0e", dash="dash"), + OptionPricingMethod.COS: ChartProps(color="#2ca02c", dash="dot"), + }, + description="Chart properties for each pricing method", + ) + + def _implied_vols( + self, r: OptionPricingResult, log_strikes: np.ndarray, ttm: float + ) -> np.ndarray: + call = np.asarray(r.call_at(log_strikes)) + intrinsic = np.maximum(0.0, 1.0 - np.exp(log_strikes)) + call = np.clip(call, intrinsic, 1.0) + return implied_black_volatility(log_strikes, call, ttm, 0.5, 1.0).values + + def _iv_error( + self, + r: OptionPricingResult, + ref: OptionPricingResult, + log_strikes: np.ndarray, + ttm: float, + ) -> float: + iv = implied_black_volatility( + log_strikes, np.asarray(r.call_at(log_strikes)), ttm, 0.5, 1.0 + ).values + iv_ref = implied_black_volatility( + log_strikes, np.asarray(ref.call_at(log_strikes)), ttm, 0.5, 1.0 + ).values + finite = np.isfinite(iv) & np.isfinite(iv_ref) + return float(np.max(np.abs(iv[finite] - iv_ref[finite]))) + + def run_ttm(self) -> None: + for ttm in self.ttms: + ms = self.model.marginal(ttm) + max_log_strike = self.max_moneyness * np.sqrt(ttm) + log_strikes = ms.option_support( + self.ref_n + 1, max_log_strike=max_log_strike + ) + ref = ms.call_option(self.ref_n, max_log_strike=max_log_strike) + iv_ref = self._implied_vols(ref, log_strikes, ttm) + moneyness_ref = log_strikes / np.sqrt(ttm) + ttm_label = f"TTM={ttm}" + slug = ttm_label.lower().replace("=", "").replace(".", "_") + + fig = make_subplots( + rows=1, + cols=2, + subplot_titles=( + f"Implied vol smile (N={self.ref_n})", + f"Max implied vol error vs N, reference = Lewis N={self.ref_n}", + ), + ) + fig.add_trace( + go.Scatter( + x=moneyness_ref, + y=iv_ref, + name=f"reference (N={self.ref_n})", + mode="lines", + line=dict(color="white", width=2), + legendgroup="reference", + ), + row=1, + col=1, + ) + for method, props in self.charts.items(): + errors = [] + for n in self.ns: + r = ms.call_option( + n, + max_log_strike=max_log_strike, + pricing_method=method, + ) + errors.append( + min(self._iv_error(r, ref, log_strikes, ttm), self.max_iv_error) + ) + if n == self.ns[-1]: + fig.add_trace( + go.Scatter( + x=moneyness_ref, + y=self._implied_vols(r, log_strikes, ttm), + name=method.value, + mode="lines", + line=dict(color=props.color, dash=props.dash), + legendgroup=method.value, + ), + row=1, + col=1, + ) + fig.add_trace( + go.Scatter( + x=self.ns, + y=errors, + name=method.value, + mode="lines+markers", + line=dict(color=props.color, dash=props.dash), + legendgroup=method.value, + showlegend=False, + ), + row=1, + col=2, + ) + fig.update_xaxes(title_text="moneyness (k / sqrt(ttm))", row=1, col=1) + fig.update_yaxes(title_text="implied volatility", row=1, col=1) + fig.update_xaxes(title_text="N (discretization points)", row=1, col=2) + fig.update_yaxes( + title_text="max implied vol error", type="log", row=1, col=2 + ) + fig.update_layout(title=ttm_label) + fig.write_image( + assets_path(f"pricing_method_accuracy_{slug}.png"), + width=1600, + height=800, + ) + + def run_speed(self) -> None: + # --- Speed: wall-clock time per maturity call vs n --- + fig_speed = go.Figure() + ms = self.model.marginal(1.0) + for method, props in self.charts.items(): + times = [] + for n in self.ns: + t0 = time.perf_counter() + for _ in range(self.timing_reps): + ms.call_option(n, max_log_strike=1.0, pricing_method=method) + times.append((time.perf_counter() - t0) / self.timing_reps * 1000) + fig_speed.add_trace( + go.Scatter( + x=self.ns, + y=times, + name=method.value, + mode="lines+markers", + line=dict(color=props.color, dash=props.dash), + ) + ) + fig_speed.update_layout( + title=( + f"Wall-clock time per pricing call" + f" (TTM=1.0, avg over {self.timing_reps} reps)" + ), + xaxis_title="N (discretization points)", + yaxis_title="time (ms)", + ) + fig_speed.write_image(assets_path("pricing_method_speed.png")) + + # --- Accuracy vs speed trade-off --- + tradeoff_ttm = self.ttms[1] + ms = self.model.marginal(tradeoff_ttm) + max_log_strike = self.max_moneyness * np.sqrt(tradeoff_ttm) + log_strikes = ms.option_support(self.ref_n + 1, max_log_strike=max_log_strike) + ref = ms.call_option(self.ref_n, max_log_strike=max_log_strike) + fig_tradeoff = go.Figure() + for method, props in self.charts.items(): + errors, times = [], [] + for n in self.ns: + t0 = time.perf_counter() + for _ in range(self.timing_reps): + r = ms.call_option( + n, max_log_strike=max_log_strike, pricing_method=method + ) + elapsed = (time.perf_counter() - t0) / self.timing_reps * 1000 + errors.append( + min( + self._iv_error(r, ref, log_strikes, tradeoff_ttm), + self.max_iv_error, + ) + ) + times.append(elapsed) + fig_tradeoff.add_trace( + go.Scatter( + x=times, + y=errors, + name=method.value, + mode="lines+markers", + text=[str(n) for n in self.ns], + textposition="top center", + line=dict(color=props.color, dash=props.dash), + ) + ) + fig_tradeoff.update_layout( + title=f"Accuracy vs speed trade-off (TTM={tradeoff_ttm})", + xaxis_title="time (ms)", + yaxis_title="max implied vol error", + yaxis_type="log", + ) + fig_tradeoff.write_image(assets_path("pricing_method_tradeoff.png")) + + +if __name__ == "__main__": + pr = Heston.create(vol=0.5, kappa=2, sigma=0.8, rho=-0.2) + comparison = PricingMethodComparison(model=pr) + comparison.run_ttm() + comparison.run_speed() diff --git a/docs/glossary.md b/docs/glossary.md index 28a17454..f3d7aed8 100644 --- a/docs/glossary.md +++ b/docs/glossary.md @@ -25,6 +25,14 @@ The characteristic exponent $\phi_{x,u}$ is defined from the \Phi_{x,u} = e^{-\phi_{x,u}} \end{equation} +The library implements the [characteristic_exponent][quantflow.sp.base.StochasticProcess1D.characteristic_exponent] for several stochastic processes, +including Brownian motion, Poisson and compound Poisson processes, the CIR square-root +diffusion, Ornstein-Uhlenbeck processes, Heston and Double Heston stochastic volatility +models, jump-diffusion models, and the Barndorff-Nielsen-Shephard (BNS) model. +Having an analytic form of the characteristic exponent for these processes +enables efficient option pricing via Fourier inversion methods such as the +[Lewis (2001)](bibliography.md#lewis) and [Carr-Madan (1999)](bibliography.md#carr_madan) approaches. + ## Cumulative Distribution Function (CDF) The cumulative distribution function (CDF), or just distribution function, @@ -81,7 +89,7 @@ The log-strike is used as input for all Black-Scholes type formulas. Moneyness is used in the context of option pricing in order to compare options with different maturities. It is defined as \begin{equation} - m = \frac{1}{\sqrt{\tau}}\ln{\frac{K}{F}} + m = \frac{1}{\sqrt{\tau}}\ln{\frac{K}{F}} = \frac{k}{\sqrt{\tau}} \end{equation} where $K$ is the strike, $F$ is the Forward price, and $\tau$ is the time to maturity. It is used to compare options with different maturities by scaling the [log-strike](#log-strike) by the square root of time to maturity. This is because the price of the underlying asset is subject to random fluctuations, if these fluctuations follow a Brownian motion than the standard deviation of the price movement will increase with the square root of time. @@ -97,6 +105,24 @@ The vol-adjusted moneyness is used in the context of option pricing in order to where $K$ is the strike, $F$ is the Forward price, $\tau$ is the time to maturity and $\sigma$ is the implied Black volatility. +## Parseval's Theorem + +Parseval's theorem states that for two square-integrable functions $f$ and $g$ with Fourier transforms $\hat{f}$ and $\hat{g}$ + +\begin{equation} +\int_{-\infty}^\infty f(s)\, g(s)\, ds = \frac{1}{2\pi} \int_{-\infty}^\infty \hat{f}(u)\, \overline{\hat{g}(u)}\, du +\end{equation} + +where $\overline{\hat{g}(u)}$ denotes the complex conjugate of $\hat{g}(u)$. + +If $g$ is real-valued, then + +\begin{equation} +\overline{\hat{g}(u)} = \hat{g}(-u) +\end{equation} + +It is used in the derivation of the [Lewis option pricing formula](theory/option_pricing.md#lewis-formula). + ## Probability Density Function (PDF) The [probability density function](https://en.wikipedia.org/wiki/Probability_density_function) diff --git a/docs/javascripts/mathjax.js b/docs/javascripts/mathjax.js index 80e81ba5..6bc413c2 100644 --- a/docs/javascripts/mathjax.js +++ b/docs/javascripts/mathjax.js @@ -2,6 +2,7 @@ window.MathJax = { tex: { inlineMath: [["\\(", "\\)"]], displayMath: [["\\[", "\\]"]], + tags: "ams", processEscapes: true, processEnvironments: true }, diff --git a/docs/references.bib b/docs/references.bib index f414dbbc..877008dd 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -2,138 +2,151 @@ --- @book{bertoin, - title={Lévy Processes}, - author={Bertoin, J.}, - year={1996}, - publisher={Cambridge University Press}, + title={Lévy Processes}, + author={Bertoin, J.}, + year={1996}, + publisher={Cambridge University Press}, } -@article{carr_madan, - title = {Option valuation using the fast Fourier transform}, - journal = {Journal of Computational Finance}, - volume = {3}, - pages = {463-520}, - author = {Carr, P. and Madan, D.}, - year = {1999}, - url = {http://faculty.baruch.cuny.edu/lwu/890/CarrMadan99.pdf}, +@article{bns, + title={Non-Gaussian OU based models and some of their uses in financial economics}, + author={Barndorff-Nielsen, O.E. & Shephard, N.}, + journal={Journal of the Royal Statistical Society}, + year={2001}, + volume={63}, + number={2}, + url={https://drive.google.com/file/d/11TnGl2ER_-QzEu_h7ZtGGass2rF8_u1Y/view}, } -@article{lee_option, - title = {Option Pricing by Transform Methods: Extensions, Unification, and Error Control}, - author = {Lee, R. W.}, - year = {2004}, - url = {https://www.math.uchicago.edu/~rl/dft.pdf}, +@article{carr_madan, + title={Option valuation using the fast Fourier transform}, + journal={Journal of Computational Finance}, + volume={3}, + pages={463-520}, + author={Carr, P. and Madan, D.}, + year={1999}, + url={http://faculty.baruch.cuny.edu/lwu/890/CarrMadan99.pdf}, } @article{carr_wu, - title = {Time-changed Lévy processes and option pricing}, - author = {Carr, P. and Wu, L.}, - year = {2002}, - journal = {Journal of Financial Economics}, - volume = {7}, - pages = {113-141}, - url = {https://engineering.nyu.edu/sites/default/files/2019-03/Carr-time-changed-levy-processes-option-pricing.pdf}, + title={Time-changed Lévy processes and option pricing}, + author={Carr, P. and Wu, L.}, + year={2002}, + journal={Journal of Financial Economics}, + volume={7}, + pages={113-141}, + url={https://engineering.nyu.edu/sites/default/files/2019-03/Carr-time-changed-levy-processes-option-pricing.pdf}, } @article{cgmy, - title={Stochastic Volatility for Lévy processes}, - author={Carr P., Geman H., Madan D.B. and Yor M.}, - year={2003}, - journal = {Mathematical Finance}, - volume = {13}, - number = {3}, - url = {https://engineering.nyu.edu/sites/default/files/2019-03/Carr-stochastic-volatility-levy-processes.pdf}, + title={Stochastic Volatility for Lévy processes}, + author={Carr P., Geman H., Madan D.B. and Yor M.}, + year={2003}, + journal={Mathematical Finance}, + volume={13}, + number={3}, + url={https://engineering.nyu.edu/sites/default/files/2019-03/Carr-stochastic-volatility-levy-processes.pdf}, } -@article{frft, +@article{chourdakis, title={Option Pricing Using the Fractional FFT}, author={Chourdakis, K.}, journal={Journal of Computational Finance}, year={2004}, volume={8}, pages={1-18}, - url = {https://citeseerx.ist.psu.edu/document?repid=rep1&type=pdf&doi=6bdf4696312d37427eda2740137650c09deacda7} + url={https://citeseerx.ist.psu.edu/document?repid=rep1&type=pdf&doi=6bdf4696312d37427eda2740137650c09deacda7}, } -@mastersthesis{saez, - title={Fourier Transform Methods for Option Pricing: An Application to extended Heston-type Models}, - author={Saez, G. K. G.}, - school={Universidad del Pais Vasco}, - year={2014}, - url = {https://www.uv.es/bfc/TFM2014/008-014.pdf} +@article{cos, + title={A Novel Pricing Method for European Options Based on Fourier-Cosine Series Expansions}, + author={Fang, F. and Oosterlee, K.}, + year={2008}, + url={https://mpra.ub.uni-muenchen.de/8914/4/MPRA_paper_8914.pdf}, } -@article{heston-calibration, - url = {https://doi.org/10.1515/math-2017-0058}, - title = {Calibration and simulation of Heston model}, - author = {Milan Mrázek and Jan Pospíšil}, - pages = {679--704}, - volume = {15}, - number = {1}, - journal = {Open Mathematics}, - doi = {doi:10.1515/math-2017-0058}, - year = {2017}, - lastchecked = {2023-07-07} +@article{dsp, + title={Doubly Stochastic Poisson Processes with Affine Intensities}, + author={Unknown}, + journal={internet}, + year={2017}, + url={https://editorialexpress.com/cgi-bin/conference/download.cgi?db_name=sbe35&paper_id=179}, } -@article{heston-simulation, - url = {https://papers.ssrn.com/sol3/papers.cfm?abstract_id=946405}, - title = {Efficient Simulation of the Heston Stochastic Volatility Model}, - author = {Andersen, Leif B.G.}, - volume = {11}, - number = {3}, - journal = {Journal of Computational Finance}, - year = {2008}, +@article{ekf, + title={Parameter Estimations of Heston Model Based on Consistent Extended Kalman Filter}, + author={Wang X., He X., Zhao Y. & Zuo Z.}, + journal={internet}, + year={2017}, + url={https://www.sciencedirect.com/science/article/pii/S2405896317324758}, } -@article{ukf, - title={The Unscented Kalman Filter for Nonlinear Estimation}, - author={Merwe}, - journal={internet}, - year={2014}, - url={https://groups.seas.harvard.edu/courses/cs281/papers/unscented.pdf} +@article{gamma-ou, + title={Gamma Related Ornstein-Uhlenbeck Processes and their Simulation}, + author={Sabino, P. & Petroni, C.}, + journal={Journal of Statistical Computation and Simulation}, + year={2021}, + volume={91}, + number={6}, + url={https://doi.org/10.1080/00949655.2020.1842408}, } -@article{ekf, - title={Parameter Estimations of Heston Model Based on Consistent Extended Kalman Filter}, - author={Wang X., He X., Zhao Y. & Zuo Z.}, - journal={internet}, - year={2017}, - url={https://www.sciencedirect.com/science/article/pii/S2405896317324758}, +@article{heston-calibration, + title={Calibration and simulation of Heston model}, + author={Milan Mrázek and Jan Pospíšil}, + journal={Open Mathematics}, + year={2017}, + volume={15}, + number={1}, + pages={679--704}, + doi={doi:10.1515/math-2017-0058}, + url={https://doi.org/10.1515/math-2017-0058}, } -@article{ou, - title={Non-Gaussian OU based models and some of their uses in financial economics}, - author={Barndorff-Nielsen, O.E. & Shephard, N.}, - journal={Journal of the Royal Statistical Society}, - year={2001}, - volume = {63}, - number = {2}, - url={https://www.sciencedirect.com/science/article/pii/S2405896317324758}, +@article{heston-simulation, + title={Efficient Simulation of the Heston Stochastic Volatility Model}, + author={Andersen, Leif B.G.}, + journal={Journal of Computational Finance}, + year={2008}, + volume={11}, + number={3}, + url={https://papers.ssrn.com/sol3/papers.cfm?abstract_id=946405}, } -@article{gamma-ou, - title="Gamma Related Ornstein-Uhlenbeck Processes and their Simulation", - author={Sabino, P. & Petroni, C.}, - journal={Journal of Statistical Computation and Simulation}, - year={2021}, - volume = {91}, - number = {6}, - url={https://doi.org/10.1080/00949655.2020.1842408}, +@article{lee_option, + title={Option Pricing by Transform Methods: Extensions, Unification, and Error Control}, + author={Lee, R. W.}, + year={2004}, + url={https://www.math.uchicago.edu/~rl/dft.pdf}, } -@article{dspp, - url={https://editorialexpress.com/cgi-bin/conference/download.cgi?db_name=sbe35&paper_id=179}, - title={Doubly Stochastic Poisson Processes with Affine Intensities}, - author={Unknown}, - journal={internet}, - year={2017}, +@article{lewis, + title={A Simple Option Formula for General Jump-Diffusion and other Exponential Lévy Processes}, + author={Lewis, A. L.}, + year={2001}, + url={https://drive.google.com/file/d/1JiYfyOU7lUKHrMqTRPbBCdasctwYDc4Z/view}, } @mastersthesis{molnar, - url={https://drive.google.com/file/d/1zCU1OZyrKQLpxaypPv9U5UPbReBDXcMf/view}, - title={Volatility modeling and forecasting: utilization of realized volatility, implied volatility and the highest and lowest price of the day}, - author={Peter Molnar}, - school={University of Economics in Prague}, - year={2020}, + title={Volatility modeling and forecasting: utilization of realized volatility, implied volatility and the highest and lowest price of the day}, + author={Peter Molnar}, + school={University of Economics in Prague}, + year={2020}, + url={https://drive.google.com/file/d/1zCU1OZyrKQLpxaypPv9U5UPbReBDXcMf/view}, +} + +@mastersthesis{saez, + title={Fourier Transform Methods for Option Pricing: An Application to extended Heston-type Models}, + author={Saez, G. K. G.}, + school={Universidad del Pais Vasco}, + year={2014}, + url={https://www.uv.es/bfc/TFM2014/008-014.pdf}, +} + +@article{ukf, + title={The Unscented Kalman Filter for Nonlinear Estimation}, + author={Merwe}, + journal={internet}, + year={2014}, + url={https://groups.seas.harvard.edu/courses/cs281/papers/unscented.pdf}, } diff --git a/docs/theory/characteristic.md b/docs/theory/characteristic.md index 94ae1689..73184aa0 100644 --- a/docs/theory/characteristic.md +++ b/docs/theory/characteristic.md @@ -11,6 +11,12 @@ The characteristic function of a random variable $x$ is the Fourier (inverse) tr \Phi_{x,u} = {\mathbb E}\left[e^{i u x}\right] = \int e^{i u s} {\mathbb P}_x\left(ds\right) \end{equation} +If the distribution of $x$ has a [probability density](../glossary.md#probability-density-function-pdf) $f_x$, then the characteristic function is the Fourier transform of $f_x$: + +\begin{equation} + \Phi_{x,u} = \int e^{i u s} f_x\left(s\right) ds +\end{equation} + ## Properties * $\Phi_{x, 0} = 1$ diff --git a/docs/theory/convexity_correction.md b/docs/theory/convexity_correction.md new file mode 100644 index 00000000..4eb4a7ea --- /dev/null +++ b/docs/theory/convexity_correction.md @@ -0,0 +1,119 @@ +# Convexity Correction + +When pricing derivatives we work under a risk-neutral measure and require that the +discounted asset price is a martingale. If we model the log-return directly as a +stochastic process $x_t$, that process will generally not satisfy this condition on +its own. The convexity correction is the deterministic adjustment that restores it. + +## Setup + +Consider an asset price of the form + +\begin{equation} +S_t = S_0\, e^{s_t} +\end{equation} + +where $s_t$ is the log-return. We model $s_t$ as + +\begin{equation} +s_t = x_t - c_t +\end{equation} + +where $x_t$ is the driving stochastic process (e.g. a Brownian motion, a Lévy process, +or a stochastic-volatility process) and $c_t$ is a deterministic function of time. + +The no-arbitrage condition (with zero interest rates) requires the forward price to +equal the spot price, i.e. ${\mathbb E}[S_t] = S_0$, which is equivalent to + +\begin{equation} +{\mathbb E}\!\left[e^{s_t}\right] = 1 +\end{equation} + +## The Correction Term + +Substituting $s_t = x_t - c_t$ into the martingale condition gives + +\begin{align*} +{\mathbb E}\!\left[e^{x_t - c_t}\right] &= 1 \\ +e^{-c_t}\,{\mathbb E}\!\left[e^{x_t}\right] &= 1 +\end{align*} + +so the correction must satisfy + +\begin{equation} +c_t = \log {\mathbb E}\!\left[e^{x_t}\right] +\end{equation} + +In terms of the [characteristic exponent](characteristic.md) $\phi_{x_t, u}$ (defined by +$\Phi_{x_t}(u) = e^{-\phi_{x_t,u}}$), evaluating at $u = -i$ gives + +\begin{equation} +{\mathbb E}\!\left[e^{x_t}\right] = \Phi_{x_t}(-i) = e^{-\phi_{x_t,-i}} +\end{equation} + +therefore + +\begin{equation} +c_t = -\phi_{x_t, -i} +\end{equation} + +This is the log-Laplace transform (cumulant generating function) of $x_t$ evaluated +at 1. It is always real-valued and non-negative by Jensen's inequality, since +$e^{x_t}$ is convex (hence the name *convexity correction*). + +## Effect on the Characteristic Function + +The characteristic function of the log-return $s_t = x_t - c_t$ is + +\begin{equation} +\Phi_{s_t}(u) = {\mathbb E}\!\left[e^{iu s_t}\right] = {\mathbb E}\!\left[e^{iu(x_t - c_t)}\right] = e^{-iuc_t}\,\Phi_{x_t}(u) +\end{equation} + +The correction introduces a phase shift proportional to $c_t$. Two key values confirm +the martingale property: + +* $\Phi_{s_t}(0) = 1$ (normalization, always true) +* $\Phi_{s_t}(-i) = e^{-c_t}\,{\mathbb E}[e^{x_t}] = e^{-c_t} e^{c_t} = 1$ (martingale condition) + +## Wiener Process + +For a Brownian motion $x_t = \sigma W_t$ the moment generating function is + +\begin{equation} +{\mathbb E}\!\left[e^{x_t}\right] = e^{\sigma^2 t / 2} +\end{equation} + +so the correction is + +\begin{equation} +c_t = \frac{\sigma^2 t}{2} +\end{equation} + +This is the classical term that appears in the Black-Scholes formula. It grows linearly +with time, reflecting the fact that geometric Brownian motion drifts upward on average +(the asset price is log-normally distributed with positive variance). + +```python +from quantflow.sp.weiner import WeinerProcess +pr = WeinerProcess(sigma=0.5) +-pr.characteristic_exponent(1, complex(0,-1)) # c_t at t=1 +``` + +which is the same as + +```python +pr.convexity_correction(1) +``` + +## General Lévy Process + +For any Lévy process $x_t$ with characteristic exponent $\psi$ (so that +$\phi_{x_t, u} = t\,\psi(u)$), the correction is + +\begin{equation} +c_t = -t\,\psi(-i) = t\,\log {\mathbb E}\!\left[e^{x_1}\right] +\end{equation} + +The correction is linear in time for all Lévy processes. For processes with jumps +(e.g. Merton jump-diffusion, variance gamma) the jump component contributes an +additional positive term on top of the diffusion $\sigma^2 t / 2$. diff --git a/docs/theory/inversion.md b/docs/theory/inversion.md index fb1c37e9..8849aab3 100644 --- a/docs/theory/inversion.md +++ b/docs/theory/inversion.md @@ -21,11 +21,11 @@ f(x) &\approx \frac{1}{\pi}\sum_{m=0}^{N-1} h_m e^{-i u_m x} \Phi_x\left(u_m\rig * $N$ is the number of discretization points and must be large enough so that the characteristic function is virtually 0 for $u_{N-1}=\delta_u N$. * $h_m$ is given by the integration methodology, either trapezoidal or Simpson rule (the library supports both, with trapezoidal as default). -For full details, see Carr and Madan (1999) and Chourdakis (2005). +For full details, see [Carr & Madan](../bibliography.md#carr_madan) and [Chourdakis](../bibliography.md#chourdakis). One could use the inverse Fourier transform to solve the integral. However, it has $O(N^2)$ time complexity. One alternative, implemented in the library, is using the Fast Fourier Transform (FFT), which has $O(N \log N)$ time complexity. -Another more flexible alternative is the Fractional FFT (FRFT), described in Chourdakis (2005). This is the methodology used by default in the library. +Another more flexible alternative is the Fractional FFT (FRFT), described in [Chourdakis](../bibliography.md#chourdakis). This is the methodology used by default in the library. ## FFT Integration diff --git a/docs/theory/option_pricing.md b/docs/theory/option_pricing.md index 87d86ab6..fb6884c0 100644 --- a/docs/theory/option_pricing.md +++ b/docs/theory/option_pricing.md @@ -1,114 +1,232 @@ # Option Pricing -We can use the tooling from characteristic function inversion to price European call options on an underlying $S_t = S_0 e^{s_t}$, where $S_0$ is the spot price at time 0. +We use characteristic function inversion to price European call options on an underlying $S_t = S_0 e^{s_t}$, where $S_0$ is the spot price at time 0. We assume zero interest rates, so the forward equals the spot. The log-return $s_t = x_t - c_t$ is constructed from a driving process $x_t$ and a deterministic [convexity correction](convexity_correction.md) $c_t$ that enforces the martingale condition ${\mathbb E}[e^{s_t}] = 1$. -## Convexity Correction +## Call Option + +The price $C$ of a call option with strike $K$ is defined as + +\begin{equation} +\begin{aligned} +C &= S_0 c_k \\ +k &= \ln\frac{K}{S_0} \\ +c_k &= {\mathbb E}\left[\left(e^{s_t} - e^k\right)^+\right] + = \int_{-\infty}^\infty \left(e^s - e^k\right)^+ f_{s_t}(s)\, ds +\end{aligned} +\label{call-price} +\end{equation} + +$k$ is the [log-strike](../glossary.md#log-strike) and $f_{s_t}$ is the probability density function of $s_t$. The call price is the discounted expected payoff under the risk-neutral measure, which simplifies to the undiscounted expected payoff when interest rates are zero. -We assume an interest rate of 0, so that the forward price equals the spot price. This assumption leads to the following no-arbitrage condition +All three methods share this starting point. They all express $c_k$ via the characteristic function $\Phi_{s_t}$, but differ in how the integration contour is chosen, how the payoff is handled, and the discretisation strategy. -\begin{align} -s_t &= x_t - c_t \\ -{\mathbb E}_0\left[e^{s_t} \right] &= {\mathbb E}_0\left[e^{x_t - c_t} \right] = e^{-c_t} {\mathbb E}_0\left[e^{x_t} \right] = e^{-c_t} e^{-\phi_{x_t, -i}} = 1 -\end{align} +## Carr & Madan -Therefore, $c_t$ represents the so-called convexity correction term, and it is equal to +We follow [Carr & Madan](../bibliography.md#carr_madan) and write the Fourier transform of the call option as \begin{equation} - c_t = -\phi_{x_t, -i} +\Psi_u = \int_{-\infty}^\infty e^{i u k} c_k dk \end{equation} -The characteristic function of $s_t$ is given by +Note that $c_k$ tends to $e^{s_t}$ as $k \to -\infty$, therefore the call price function is not square-integrable. In order to obtain integrability, Carr-Madan introduces a damping factor $e^{\alpha k}$ and works with the modified call $\tilde{c}_k = e^{\alpha k} c_k$. The free parameter $\alpha > 0$ is chosen so that $\tilde{c}_k$ is square-integrable. Taking the Fourier transform of $\tilde{c}_k$ and evaluating at $u = v - i\alpha$ (real $v$) gives \begin{equation} - \Phi_{s_t}\left(u\right) = \Phi_x\left(u\right) e^{-i u c_t} +\Psi_u = \frac{\Phi_{s_t}\left(u-i\right)}{iu \left(1 + iu\right)} \end{equation} -As you can see the convexity correction increases with time horizon. Let us take a few examples: +Inverting, the call price is recovered as -### Weiner Process +\begin{equation} +\begin{aligned} +c_k &= \int_0^{\infty} e^{-iuk} \Psi\left(u\right) du \\ + &= \frac{e^{-\alpha k}}{\pi} \int_0^{\infty} e^{-ivk} \Psi\left(v-i\alpha\right) dv \\ +\end{aligned} +\end{equation} + +The same FFT/FRFT machinery used for PDF inversion applies here. + +### Choice of $\alpha$ -This is the well-known convexity correction which appears in all diffusion-driven SDEs: +For $\Psi_u$ to be well-defined, the characteristic function $\Phi_{s_t}(u - i)$ must be finite at $u = v - i\alpha$, i.e. $\Phi_{s_t}(v - i(\alpha + 1))$ must be finite. Denoting the upper edge of the strip of analyticity by $\beta_+$, this requires \begin{equation} - c_t = \frac{\sigma^2 t}{2} +\alpha + 1 < \beta_+ \end{equation} -```python -from quantflow.sp.weiner import WeinerProcess -pr = WeinerProcess(sigma=0.5) --pr.characteristic_exponent(1, complex(0,-1)) -``` +Positive values of $\alpha$ assist the integrability of the modified call value over the +negative moneyness axis, but aggravate the same condition for the positive moneyness axis. For the modified call value to be integrable in the positive moneyness +direction, and hence for it to be square-integrable as well, a sufficient condition +is provided by $\Psi_{-i\alpha}$ being finite, which means the characteristic function $\Phi_{t,{-(\alpha+1)i}}$ is finite. -which is the same as +The constraint $\alpha + 1 < \beta_+$ can be restrictive: if the model has thin tails then $\beta_+$ is close to 1 and no positive $\alpha$ satisfies it cleanly. A poor choice of $\alpha$ leads to numerical instabilities, especially at short maturities. -```python -pr.convexity_correction(1) -``` +## Lewis Formula -## Call Option +[Lewis](../bibliography.md#lewis) starts from the same call price expression but avoids the damping parameter by applying [Parseval's theorem](../glossary.md#parsevals-theorem) directly to $\eqref{call-price}$ and +shifting the integration in the complex plane to ensure the call options transform is well-defined and integrable. +The shift uses the [Residue theorem](https://en.wikipedia.org/wiki/Residue_theorem) to account for the poles of the payoff transform. -The price $C$ of a call option with strike $K$ is defined as -\begin{align} -C &= S_0 c_k \\ -k &= \ln\frac{K}{S_0}\\ -c_k &= {\mathbb E}\left[\left(e^{s_t} - e^k\right)1_{s_t\ge k}\right] -\end{align} +### Fourier Transform of the Payoff -We follow Carr and Madan (1999) and write the Fourier transform of the call option as +The Fourier transform of the call payoff $g(s) = (e^s - e^k)^+$ is \begin{equation} -\Psi_u = \int_{-\infty}^\infty e^{i u k} c_k dk +\hat{g}(u) = \int_{-\infty}^\infty e^{ius} (e^s - e^k)^+\, ds = \int_k^\infty e^{ius}(e^s - e^k)\, ds \end{equation} -Note that $c_k$ tends to $e^{s_t}$ as $k \to -\infty$, therefore the call price function is not square-integrable. In order to obtain integrability, we choose complex values of $u$ of the form +For the first term to be integrable, we need $\mathrm{Im}(u) > 1$, while +for the second term we need $\mathrm{Im}(u) > 0$. Therefore the Fourier +transform of the payoff is well-defined in the strip $\mathrm{Im}(u) > 1$. + +Splitting and integrating term by term gives + \begin{equation} -u = v - i \alpha +\hat{g}(u) = \frac{e^{(1+iu)k}}{iu(1+iu)} \qquad \mathrm{Im}(u) > 1 +\label{payoff-ft} \end{equation} -The value of $\alpha$ is a numerical choice we can check later. -It is possible to obtain the analytical expression of $\Psi_u$ in terms of the characteristic function $\Phi_s$. Once we have that expression, we can use the Fourier transform tooling presented previously to calculate option prices in this way +The denominator is zero at $u = 0$ and $u = i$. -\begin{align} -c_k &= \int_0^{\infty} e^{-iuk} \Psi\left(u\right) du \\ - &= \frac{e^{-\alpha k}}{\pi} \int_0^{\infty} e^{-ivk} \Psi\left(v-i\alpha\right) dv \\ -\end{align} +This result is typical, option payoffs have a Fourier transform as long as we admit a complex valued transform variable and integrate along a horizontal line in the complex plane. -The analytical expression of $\Psi_u$ is given by +### Derivation + +The call option can now be expressed via Parseval's theorem on a horizontal line in the complex plane. +We therefore evaluate the pricing integral along a horizontal contour \begin{equation} -\Psi_u = \frac{\Phi_{s_t}\left(u-i\right)}{iu \left(iu + 1\right)} + u = v + i\eta \qquad v,\eta \in \mathbb{R} \end{equation} -To integrate, we use the same approach as the PDF integration. +Since $\Phi$ is real-valued, we can replace $\overline{\Phi(u)}$ by $\Phi(-u)$: -### Choice of $\alpha$ +\begin{equation} +\begin{aligned} +c_k &= \int_{-\infty}^\infty g(s) f(s) ds \\ +&= \frac{1}{2\pi} \int_{-\infty}^\infty \hat{g}(u) \Phi_{s_t}(-u) du \\ +&= \frac{e^k}{2\pi} \int_{-\infty}^\infty \frac{e^{iuk}}{iu(1+iu)} \Phi_{s_t}(-u) du +\end{aligned} +\end{equation} -Positive values of $\alpha$ assist the integrability of the modified call value over the -negative moneyness axis, but aggravate the same condition for the positive moneyness axis. For the modified call value to be integrable in the positive moneyness -direction, and hence for it to be square-integrable as well, a sufficient condition -is provided by $\Psi_{-i\alpha}$ being finite, which means the characteristic function $\Phi_{t,{-(\alpha+1)i}}$ is finite. +The integral can be evaluated for any $\eta > 1$ as discussed in the previous section. +However, the integrand can now be computed for smaller $\eta$, provided the two `poles` at $u = 0$ and $u = i$ are avoided. -## Black Formula +If we move the contour in $\eta \in (0,1)$, the integrand is well-defined. By the residue theorem, the value of +the integral is changed by the contribution of minus $2\pi i$ times the residues at $u = i$. -Here we illustrate how to use the characteristic function integration with the classical [Weiner process](https://en.wikipedia.org/wiki/Wiener_process). +\begin{equation} +c_k = \frac{e^k}{2\pi} \int_{-\infty}^\infty \frac{e^{iuk}}{iu(1+iu)} \Phi_{s_t}(-u) du -2\pi i\, \text{Res}(i) \qquad \eta \in (0,1) +\end{equation} -```python -from quantflow.sp.weiner import WeinerProcess -from quantflow.options.bs import black_call +The residue at $u = i$ is $\frac{i}{2\pi}$, so the call price can be expressed as -ttm = 1 -p = WeinerProcess(sigma=0.5) -m = p.marginal(ttm) -m.std() -``` +\begin{equation} +c_k = 1 + \frac{e^k}{2\pi} \int_{-\infty}^\infty \frac{e^{iuk}}{iu(1+iu)} \Phi_{s_t}(-u) du \qquad \eta \in (0,1) +\end{equation} + +By choosing $\eta=\frac{1}{2}$, symmetrically located between the two poles, we can avoid numerical instabilities and obtain a stable formula for the call price: + +\begin{equation} +c_k = 1 + \frac{e^k}{2\pi} \int_{-\infty}^\infty \frac{e^{iuk}}{iu(1+iu)} \Phi_{s_t}(-u) du \qquad u = v + \frac{i}{2} +\end{equation} + +An substituting gives the final formula: + +\begin{equation} +c_k = 1 + \frac{e^k}{2\pi} \int_{0}^\infty \frac{e^{i(v + i/2)k}}{i(v + i/2)(1+i(v + i/2))} \Phi_{s_t}(-v - \frac{i}{2}) dv +\end{equation} + + +### Comparison with Carr & Madan + +The plots below use a Heston model to compare the two formulas and highlight where the +choice of $\alpha$ matters. + +At a normal maturity (TTM=0.5) both methods agree closely with auto-selected $\alpha$: + +![Carr-Madan vs Lewis prices](../assets/examples/lewis_vs_madan_prices.png) + +At short maturities the auto-selected $\alpha$ in Carr & Madan can produce significant +errors, while Lewis remains stable: + +![Short maturity difference](../assets/examples/lewis_vs_madan_short_ttm.png) + +The sensitivity to the choice of $\alpha$ is most visible at TTM=0.02. A poor choice +(e.g. $\alpha=0.25$) yields completely wrong prices deep OTM, while Lewis (dotted) is +the stable reference: + +![Alpha sensitivity](../assets/examples/lewis_vs_madan_alpha.png) + +The example code that generates these plots: ```python -import plotly.express as px -import plotly.graph_objects as go - -r = m.call_option(64, max_moneyness=1.0) -b = black_call(r.x, p.sigma, ttm) -fig = px.line(x=r.x, y=r.y, markers=True, labels=dict(x="moneyness", y="call price")) -fig.add_trace(go.Scatter(x=r.x, y=b, name="analytical", line=dict())) -fig.show() +--8<-- "docs/examples/carr_madan_vs_lewis.py" ``` + +## COS Method + +The [COS method](../bibliography.md#cos) (Fang and Oosterlee, 2008) takes a different +approach: rather than expressing the call price as a contour integral of the +characteristic function, it approximates the density $f_{s_t}$ on a truncated interval +$[a, b]$ using a cosine series and evaluates the payoff integral analytically against +each basis function. + +### Density approximation + +On the truncated interval $[a, b]$, the density is expanded as + +\begin{equation} + f_{s_t}(x) \approx \frac{2}{b-a} + \sum_{j=0}^{N-1}{}^{\prime} + \mathrm{Re}\!\left[ + \hat{\Phi}\!\left(\frac{j\pi}{b-a}\right) e^{-ij\pi a/(b-a)} + \right] + \cos\!\left(\frac{j\pi(x-a)}{b-a}\right) +\end{equation} + +where $\hat{\Phi}$ is the martingale-corrected characteristic function and the prime +denotes a half weight on the $j=0$ term. The truncation interval is chosen as +$[a, b] = [-L\sigma, L\sigma]$ where $\sigma$ is the standard deviation of $s_t$ and +$L$ is a parameter (default 12). + +### Payoff coefficients + +Substituting the density approximation into the call price integral and changing +variables to $y = x - k$ turns the payoff into the strike-independent form +$(e^y - 1)^+$. The payoff integral against each cosine basis function evaluates +analytically as + +\begin{equation} + V_j = \frac{2}{b-a}\left[\chi_j(0,b) - \psi_j(0,b)\right] +\end{equation} + +where $\nu_j = j\pi/(b-a)$ and + +\begin{equation} + \chi_j(0,b) = \int_0^b e^y \cos\!\left(\nu_j(y-a)\right) dy + = \frac{(-1)^j e^b - \cos(\nu_j a) + \nu_j\sin(\nu_j a)}{1 + \nu_j^2} +\end{equation} + +\begin{equation} + \psi_j(0,b) = \int_0^b \cos\!\left(\nu_j(y-a)\right) dy + = \begin{cases} b & j = 0 \\ \dfrac{\sin(\nu_j a)}{\nu_j} & j > 0 \end{cases} +\end{equation} + +### Call price + +Combining the density approximation with the payoff coefficients and accounting for the +change of variables gives the call price in forward space as + +\begin{equation} + c_k \approx e^k \sum_{j=0}^{N-1}{}^{\prime} + \mathrm{Re}\!\left[ + \hat{\Phi}\!\left(\frac{j\pi}{b-a}\right) + e^{-ij\pi(k+a)/(b-a)} + \right] V_j +\end{equation} + +The $e^k$ prefactor converts from the strike-normalised density of $y = x - k$ back to +forward-space pricing. Unlike Carr-Madan, no damping parameter is needed. Unlike Lewis, +the sum is over the payoff domain rather than the frequency domain, which can give +faster convergence for smooth densities at the cost of a fixed truncation error +controlled by $L$. diff --git a/docs/tutorials/pricing_method_comparison.md b/docs/tutorials/pricing_method_comparison.md new file mode 100644 index 00000000..0af06ff1 --- /dev/null +++ b/docs/tutorials/pricing_method_comparison.md @@ -0,0 +1,74 @@ +# Pricing Method Comparison + +This tutorial compares the three Fourier-based pricing methods available in +[OptionPricer][quantflow.options.pricer.OptionPricer]: +[Carr-Madan](../bibliography.md#carr_madan), +[Lewis](../bibliography.md#lewis), and +[COS](../bibliography.md#cos). +All three invert the characteristic function of the log-return to price European calls, +but they differ in how they handle the payoff transform, the integration contour, and +the discretisation strategy. See [Option Pricing](../theory/option_pricing.md) for the +full derivation of all three methods. + +## Selecting the method + +Pass `method` when constructing the pricer: + +```python +from quantflow.options.pricer import OptionPricer, OptionPricingMethod +from quantflow.sp.heston import Heston + +model = Heston.create(vol=0.2) +pricer = OptionPricer(model=model, method=OptionPricingMethod.COS) +result = pricer.maturity(1.0) +``` + +## Accuracy comparison + +The charts below use a Heston model with $\sigma=0.8$, $\kappa=2$, $\rho=-0.2$, +$\text{vol}=0.5$. The reference solution is Lewis with $N=8192$. +Implied vol errors are clipped at 10% — Lewis can produce errors well above this +at low $N$, particularly at short maturities, which would otherwise dominate the scale. + +At long maturities (TTM=1.0) all three methods converge quickly: + +[![Accuracy TTM=1.0](../assets/examples/pricing_method_accuracy_ttm1_0.png)](../assets/examples/pricing_method_accuracy_ttm1_0.png){target="_blank"} + +At medium maturities (TTM=0.25) COS tends to converge faster than Carr-Madan for the +same $N$: + +[![Accuracy TTM=0.25](../assets/examples/pricing_method_accuracy_ttm0_25.png)](../assets/examples/pricing_method_accuracy_ttm0_25.png){target="_blank"} + +At very short maturities (TTM=0.02) the differences are most pronounced. Carr-Madan +can struggle with the auto-selected $\alpha$, while Lewis and COS remain stable: + +[![Accuracy TTM=0.02](../assets/examples/pricing_method_accuracy_ttm0_02.png)](../assets/examples/pricing_method_accuracy_ttm0_02.png){target="_blank"} + +## Speed comparison + +| Method | Complexity | +|---|---| +| Carr-Madan | $O(N \log N)$ via Fractional Fourier Transform | +| Lewis | $O(N \log N)$ via Fractional Fourier Transform | +| COS | $O(N^2)$ — dense $N \times N$ complex matrix-vector product | + +COS is slower because it evaluates $e^{-i j\pi(k+a)/(b-a)}$ for every combination of +strike $k$ and cosine index $j$, forming an explicit $N \times N$ matrix before +multiplying by the coefficient vector. Lewis and Carr-Madan avoid this by evaluating +the transform on a uniform grid and applying the FrFT in $O(N \log N)$. + +Wall-clock time per pricing call as a function of $N$: + +![Speed](../assets/examples/pricing_method_speed.png) + +## Accuracy vs speed trade-off + +Each point is labelled with its $N$ value (TTM=0.25): + +![Trade-off](../assets/examples/pricing_method_tradeoff.png) + +## Example code + +```python +--8<-- "docs/examples/pricing_method_comparison.py" +``` diff --git a/mkdocs.yml b/mkdocs.yml index 610fbdc1..042baf58 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -1,5 +1,5 @@ site_name: QuantFlow -site_url: https://quantmind.github.io/quantflow/ +site_url: https://quantflow.quantmind.com site_dir: app/docs repo_name: quantmind/quantflow repo_url: https://github.com/quantmind/quantflow @@ -50,6 +50,7 @@ plugins: merge_init_into_class: true docstring_section_style: spacy signature_crossrefs: true + relative_crossrefs: true show_symbol_type_heading: true show_symbol_type_toc: true extensions: @@ -76,6 +77,7 @@ nav: - Volatility Surface: api/options/vol_surface.md - Stochastic Processes: - api/sp/index.md + - BNS Process: api/sp/bns.md - CIR Process: api/sp/cir.md - Compound Poisson: api/sp/compound_poisson.md - Doubly Stochastic Poisson: api/sp/dsp.md @@ -106,10 +108,12 @@ nav: - Tutorials: - tutorials/index.md - Option Pricing: tutorials/option_pricing.md + - Pricing Method Comparison: tutorials/pricing_method_comparison.md - Volatility Surface: tutorials/volatility_surface.md - Theory: - theory/index.md - Characteristic Function: theory/characteristic.md + - Convexity Correction: theory/convexity_correction.md - Inversion: theory/inversion.md - Lévy Process: theory/levy.md - Option Pricing: theory/option_pricing.md @@ -124,6 +128,7 @@ nav: - Contributing: contributing.md - Bibliography: bibliography.md markdown_extensions: + - attr_list - tables - pymdownx.arithmatex: generic: true @@ -138,5 +143,6 @@ markdown_extensions: permalink: true extra_javascript: + - javascripts/mathjax.js - https://polyfill.io/v3/polyfill.min.js?features=es6 - https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-mml-chtml.js diff --git a/quantflow/options/calibration.py b/quantflow/options/calibration.py index 4c598570..1819d14b 100644 --- a/quantflow/options/calibration.py +++ b/quantflow/options/calibration.py @@ -33,7 +33,7 @@ class OptionEntry(BaseModel): """ ttm: float = Field(description="Time to maturity in years") - moneyness: float = Field(description="Log-moneyness: log(strike / forward)") + log_strike: float = Field(description="Log-strike as log(K/F)") options: list[OptionPrice] = Field(default_factory=list) """Bid and ask option prices for this entry""" _mid_price: float | None = PrivateAttr(default=None) @@ -112,7 +112,7 @@ def model_post_init(self, _ctx: Any) -> None: for option in self.vol_surface.option_prices(converged=True): key = ModelCalibrationEntryKey(option.maturity, option.strike) if key not in self.options: - entry = OptionEntry(ttm=option.ttm, moneyness=option.moneyness) + entry = OptionEntry(ttm=option.ttm, log_strike=option.log_strike) self.options[key] = entry entry.options.append(option) @@ -178,9 +178,9 @@ def fit(self) -> OptimizeResult: self.set_params(result.x) return result - def cost_weight(self, ttm: float, moneyness: float) -> float: - """Weight for a given time to maturity and moneyness""" - return np.exp(-self.moneyness_weight * moneyness) + def cost_weight(self, ttm: float, log_strike: float) -> float: + """Weight for a given time to maturity and log-strike""" + return np.exp(-self.moneyness_weight * log_strike) def penalize(self) -> float: """Additional scalar penalty added to the cost function (default: 0)""" @@ -193,8 +193,8 @@ def residuals(self, params: np.ndarray) -> np.ndarray: res = [] with np.errstate(all="ignore"): for entry in self.options.values(): - model_price = self.pricer.call_price(entry.ttm, entry.moneyness) - weight = self.cost_weight(entry.ttm, entry.moneyness) + model_price = self.pricer.call_price(entry.ttm, entry.log_strike) + weight = self.cost_weight(entry.ttm, entry.log_strike) r = weight * (model_price - entry.mid_price()) res.append(r if np.isfinite(r) else 1e6) return np.asarray(res) @@ -218,7 +218,7 @@ def plot( model = self.pricer.maturity(cross.ttm(self.ref_date)) if max_moneyness_ttm is not None: model = model.max_moneyness_ttm( - max_moneyness_ttm=max_moneyness_ttm, support=support + max_moneyness=max_moneyness_ttm, support=support ) return plot.plot_vol_surface( pd.DataFrame([d.info_dict() for d in options]), @@ -252,7 +252,7 @@ def plot_maturities( model = self.pricer.maturity(cross.ttm(self.ref_date)) if max_moneyness_ttm is not None: model = model.max_moneyness_ttm( - max_moneyness_ttm=max_moneyness_ttm, support=support + max_moneyness=max_moneyness_ttm, support=support ) plot.plot_vol_surface( pd.DataFrame([d.info_dict() for d in options]), diff --git a/quantflow/options/divfm/pricer.py b/quantflow/options/divfm/pricer.py index 0f4d5768..93fa028f 100644 --- a/quantflow/options/divfm/pricer.py +++ b/quantflow/options/divfm/pricer.py @@ -124,14 +124,14 @@ def _compute_maturity(self, ttm: float, **kwargs: Any) -> MaturityPricer: iv_grid = np.clip(F @ self.betas, 1e-6, None) - # Convert from time-scaled moneyness M to log-moneyness m = log(K/F) - moneyness = M_grid * np.sqrt(ttm) - call = np.asarray(black_call(moneyness, iv_grid, ttm=ttm)) + # Convert from time-scaled moneyness M to log-strike = log(K/F) + log_strike = M_grid * np.sqrt(ttm) + call = np.asarray(black_call(log_strike, iv_grid, ttm=ttm)) return MaturityPricer( ttm=ttm, std=float(np.mean(iv_grid) * np.sqrt(ttm)), - moneyness=moneyness, + log_strike=log_strike, call=call, name="DIVFM", ) diff --git a/quantflow/options/pricer.py b/quantflow/options/pricer.py index 30efcebd..67e8aa96 100644 --- a/quantflow/options/pricer.py +++ b/quantflow/options/pricer.py @@ -11,6 +11,7 @@ from quantflow.sp.base import StochasticProcess1D from quantflow.utils import plot +from quantflow.utils.marginal import OptionPricingMethod from quantflow.utils.numbers import DecimalNumber, to_decimal from ..utils.types import FloatArray @@ -23,8 +24,8 @@ TTM_FACTOR = 10000 -def get_intrinsic_value(moneyness: FloatArray) -> FloatArray: - return 1.0 - np.exp(np.clip(moneyness, None, 0)) +def get_intrinsic_value(log_strike: FloatArray) -> FloatArray: + return 1.0 - np.exp(np.clip(log_strike, None, 0)) class ModelOptionPrice(BaseModel, frozen=True): @@ -38,7 +39,8 @@ class ModelOptionPrice(BaseModel, frozen=True): strike: DecimalNumber = Field(description="Strike price of the option") option_type: OptionType = Field(description="Type of the option, call or put") forward: DecimalNumber = Field(description="Forward price of the underlying") - moneyness: float = Field(description="Moneyness as log(Strike/Forward)") + log_strike: float = Field(description="Log strike over forward, i.e. log(K/F)") + moneyness: float = Field(description="Moneyness") ttm: float = Field(default=0, description="Time to maturity in years") price: float = Field(description=("Price in forward space")) delta: float = Field(description="Model delta of the option") @@ -48,9 +50,9 @@ class ModelOptionPrice(BaseModel, frozen=True): @property def black(self) -> BlackSensitivities: """Calculate the Black price for the option using the price as time value and - moneyness as log-strike""" + log-strike""" return BlackSensitivities.calculate( - k=self.moneyness, + k=self.log_strike, ttm=self.ttm, s=self.option_type.call_put(), price=self.price, @@ -70,9 +72,9 @@ def intrinsic_value(self) -> float: is positive, i.e. when the strike is above the forward price. """ if self.option_type == OptionType.call: - return max(0.0, 1.0 - np.exp(self.moneyness)) + return max(0.0, 1.0 - np.exp(self.log_strike)) else: - return max(0.0, np.exp(self.moneyness) - 1.0) + return max(0.0, np.exp(self.log_strike) - 1.0) def as_option_type( self, @@ -100,16 +102,16 @@ class MaturityPricer(BaseModel, arbitrary_types_allowed=True): ttm: float = Field(description="Time to maturity") std: float = Field(description="Standard deviation model of log returns") - moneyness: FloatArray = Field(description="Moneyness as log Strike/Forward") + log_strike: FloatArray = Field(description="Log strike over forward, i.e. log(K/F)") call: FloatArray = Field( - description="Call prices in forward space for the given moneyness" + description="Call prices in forward space for the given log strike array" ) name: str = Field(default="", description="Name of the model") @property - def moneyness_ttm(self) -> FloatArray: + def moneyness(self) -> FloatArray: """Time adjusted moneyness array""" - return self.moneyness / np.sqrt(self.ttm) + return self.log_strike / np.sqrt(self.ttm) @property def time_value(self) -> FloatArray: @@ -119,7 +121,7 @@ def time_value(self) -> FloatArray: @property def intrinsic_value(self) -> FloatArray: """Intrinsic value of the option""" - return get_intrinsic_value(self.moneyness) + return get_intrinsic_value(self.log_strike) def price( self, @@ -133,26 +135,27 @@ def price( """ strike_ = to_decimal(strike) forward_ = to_decimal(forward) - moneyness = float((strike_ / forward_).ln()) + log_strike = float((strike_ / forward_).ln()) return ModelOptionPrice( option_type=OptionType.call, strike=strike_, forward=forward_, - moneyness=moneyness, + log_strike=log_strike, + moneyness=log_strike / np.sqrt(self.ttm), ttm=self.ttm, - price=self.call_price(moneyness), - delta=self.call_delta(moneyness), - gamma=self.call_gamma(moneyness), + price=self.call_price(log_strike), + delta=self.call_delta(log_strike), + gamma=self.call_gamma(log_strike), ).as_option_type(option_type) @property def implied_vols(self) -> FloatArray: """Calculate the implied volatility""" return implied_black_volatility( - self.moneyness, + self.log_strike, self.call, ttm=self.ttm, - initial_sigma=0.5 * np.ones_like(self.moneyness), + initial_sigma=0.5 * np.ones_like(self.log_strike), call_put=1.0, ).values @@ -161,53 +164,51 @@ def df(self) -> pd.DataFrame: """Get a dataframe with the results""" return pd.DataFrame( { + "log_strike": self.log_strike, "moneyness": self.moneyness, - "moneyness_ttm": self.moneyness_ttm, "call": self.call, "implied_vol": self.implied_vols, "time_value": self.time_value, } ) - def call_price(self, moneyness: float) -> float: + def call_price(self, log_strike: float) -> float: """Price a single call option""" - return float(np.interp(moneyness, self.moneyness, self.call)) + return float(np.interp(log_strike, self.log_strike, self.call)) - def call_delta(self, moneyness: float) -> float: + def call_delta(self, log_strike: float) -> float: r"""Delta of a call option as change in price per unit change in forward. Since prices are stored in forward space (c = C/F) and m = log(K/F), the chain rule gives: dC/dF = c - dc/dm """ - dc_dm = np.gradient(self.call, self.moneyness) - return float(np.interp(moneyness, self.moneyness, self.call - dc_dm)) + dc_dm = np.gradient(self.call, self.log_strike) + return float(np.interp(log_strike, self.log_strike, self.call - dc_dm)) - def call_gamma(self, moneyness: float) -> float: + def call_gamma(self, log_strike: float) -> float: """Gamma of a call option as change in delta per unit change in forward. Since prices are stored in forward space (c = C/F) and m = log(K/F), the chain rule gives: d2C/dF2 = d2c/dm2 - dc/dm """ - dc_dm = np.gradient(self.call, self.moneyness) - d2c_dm2 = np.gradient(dc_dm, self.moneyness) - return float(np.interp(moneyness, self.moneyness, d2c_dm2 - dc_dm)) + dc_dm = np.gradient(self.call, self.log_strike) + d2c_dm2 = np.gradient(dc_dm, self.log_strike) + return float(np.interp(log_strike, self.log_strike, d2c_dm2 - dc_dm)) - def interp(self, moneyness: FloatArray) -> Self: + def interp(self, log_strike: FloatArray) -> Self: """get interpolated prices""" return self.model_copy( update=dict( - moneyness=moneyness, - call=np.interp(moneyness, self.moneyness, self.call), + log_strike=log_strike, + call=np.interp(log_strike, self.log_strike, self.call), ) ) - def max_moneyness_ttm( - self, max_moneyness_ttm: float = 1.0, support: int = 51 - ) -> Self: + def max_moneyness_ttm(self, max_moneyness: float = 1.0, support: int = 51) -> Self: """Calculate the implied volatility""" - moneyness_ttm = np.linspace(-max_moneyness_ttm, max_moneyness_ttm, support) - moneyness = np.asarray(moneyness_ttm) * np.sqrt(self.ttm) - return self.interp(moneyness) + moneyness = np.linspace(-max_moneyness, max_moneyness, support) + log_strike = np.asarray(moneyness) * np.sqrt(self.ttm) + return self.interp(log_strike) def black(self) -> Self: """Calculate the Maturity Result for the Black model with same std""" @@ -215,7 +216,7 @@ def black(self) -> Self: update=dict( call=np.asarray( black_call( - self.moneyness, self.std / np.sqrt(self.ttm), ttm=self.ttm + self.log_strike, self.std / np.sqrt(self.ttm), ttm=self.ttm ) ), name="Black", @@ -230,7 +231,7 @@ def plot(self, series: str = "implied_vol", **kwargs: Any) -> Any: class OptionPricerBase(BaseModel, arbitrary_types_allowed=True): """Abstract base class for option pricers. - Provides caching of [MaturityPricer][quantflow.options.pricer.MaturityPricer] + Provides caching of [MaturityPricer][..MaturityPricer] results and generic pricing/plotting methods. Subclasses implement ``_compute_maturity`` to define how call prices are computed for a given time to maturity. @@ -284,13 +285,13 @@ def price( def call_price( self, ttm: Annotated[float, Doc("Time to maturity")], - moneyness: Annotated[float, Doc("Moneyness as log(Strike/Forward)")], + log_strike: Annotated[float, Doc("Log strike over forward, i.e. log(K/F)")], ) -> float: """Price a single call option This method will use the cache to get the maturity pricer if possible """ - return self.maturity(ttm).call_price(moneyness) + return self.maturity(ttm).call_price(log_strike) def plot3d( self, @@ -313,12 +314,12 @@ def plot3d( maturity = self.maturity(cast(float, t)) implied[i, :] = maturity.interp(moneyness_ttm * np.sqrt(t)).implied_vols properties: dict = dict( - xaxis_title="moneyness_ttm", + xaxis_title="moneyness", yaxis_title="TTM", colorscale="viridis", dragmode=dragmode, scene=dict( - xaxis=dict(title="moneyness_ttm"), + xaxis=dict(title="moneyness"), yaxis=dict(title="TTM"), zaxis=dict(title="implied_vol"), ), @@ -340,7 +341,7 @@ class OptionPricer(OptionPricerBase, Generic[M]): """Pricer for options based on a stochastic process model. Computes call prices via the inverse Fourier transform of the - characteristic function of the underlying stochastic process. + call option transform function of the underlying stochastic process. """ model: M = Field(description="Stochastic process model for the underlying") @@ -348,19 +349,25 @@ class OptionPricer(OptionPricerBase, Generic[M]): default=128, description="Number of discretization points for the marginal distribution", ) + method: OptionPricingMethod = Field( + default=OptionPricingMethod.LEWIS, + description="Method to use for option pricing", + ) max_moneyness_ttm: float = Field( - default=1.5, description="Max time-adjusted moneyness to calculate prices" + default=1.5, description="Max moneyness to calculate prices" ) def _compute_maturity(self, ttm: float, **kwargs: Any) -> MaturityPricer: marginal = self.model.marginal(ttm) + max_log_strike = self.max_moneyness_ttm * np.sqrt(ttm) transform = marginal.call_option( - self.n, max_moneyness=self.max_moneyness_ttm * np.sqrt(ttm), **kwargs + self.n, pricing_method=self.method, max_log_strike=max_log_strike, **kwargs ) + log_strike = marginal.option_support(self.n + 1, max_log_strike=max_log_strike) return MaturityPricer( ttm=ttm, std=float(np.max(marginal.std())), - moneyness=transform.x, - call=transform.y, + log_strike=log_strike, + call=transform.call_at(log_strike), name=type(self.model).__name__, ) diff --git a/quantflow/options/surface.py b/quantflow/options/surface.py index fc9447be..eca93eb9 100644 --- a/quantflow/options/surface.py +++ b/quantflow/options/surface.py @@ -315,7 +315,7 @@ def create( open_interest: Number = ZERO, volume: Number = ZERO, inverse: bool = True, - ) -> OptionPrice: + ) -> Self: """Create an option price mainly used for testing @@ -368,7 +368,8 @@ def volume(self) -> Decimal: return self.meta.volume @property - def moneyness(self) -> float: + def log_strike(self) -> float: + """Log strike of the option, calculated as log(strike/forward)""" return float(np.log(float(self.strike / self.forward))) @property @@ -433,11 +434,11 @@ def put_price(self) -> Decimal: def is_in_the_money(self, forward: Decimal) -> bool: return self.meta.is_in_the_money(forward) - def calculate_price(self) -> OptionPrice: + def calculate_price(self) -> Self: price = Decimal( sigfig( black_price( - np.asarray(self.moneyness), + np.asarray(self.log_strike), self.implied_vol, self.ttm, 1 if self.option_type.is_call() else -1, @@ -453,8 +454,8 @@ def info_dict(self) -> dict[str, Any]: strike=float(self.strike), forward=float(self.forward), maturity=self.maturity, - moneyness=self.moneyness, - moneyness_ttm=self.moneyness / np.sqrt(self.ttm), + log_strike=self.log_strike, + moneyness=self.log_strike / np.sqrt(self.ttm), ttm=self.ttm, implied_vol=self.implied_vol, price=float(self.price_in_forward_space), @@ -473,7 +474,7 @@ class OptionArrays(NamedTuple): options: list[OptionPrice] """List of option prices corresponding to the arrays below""" - moneyness: FloatArray + log_strike: FloatArray """The log strike of the options, calculated as log(strike/forward)""" price: FloatArray """The option prices""" @@ -881,7 +882,7 @@ def disable_outliers( than 20% of the mid vol. Options with a zero mid vol are also disabled. Second pass: an [SVI][quantflow.options.svi.SVI] smile is fitted to the - surviving options (mid implied vol vs log-moneyness). Options whose + surviving options (mid implied vol vs log-strike). Options whose residual from the SVI fit exceeds `svi_residual_fraction` of their mid implied vol are disabled. This is repeated up to `repeat` times, refitting after each removal. The loop stops early if no outliers are @@ -1130,7 +1131,7 @@ def bs( with warnings.catch_warnings(): warnings.simplefilter("ignore") result = implied_black_volatility( - k=d.moneyness, + k=d.log_strike, price=d.price, ttm=d.ttm, initial_sigma=d.implied_vol, @@ -1159,7 +1160,7 @@ def calc_bs_prices( otherwise the price calculation won't be correct. """ d = self.as_array(select=select, index=index, converged=True) - return black_price(k=d.moneyness, sigma=d.implied_vol, ttm=d.ttm, s=d.call_put) + return black_price(k=d.log_strike, sigma=d.implied_vol, ttm=d.ttm, s=d.call_put) def options_df( self, @@ -1209,7 +1210,7 @@ def as_array( and price calculation It returns an [OptionArrays][quantflow.options.surface.OptionArrays] instance, - which contains the option prices and their corresponding moneyness, + which contains the option prices and their corresponding log strikes, time to maturity and implied volatility in numpy arrays for efficient calculations. """ @@ -1221,20 +1222,20 @@ def as_array( converged=converged, ) ) - moneyness = [] + log_strike = [] ttm = [] price = [] vol = [] call_put = [] for option in options: - moneyness.append(float(option.moneyness)) + log_strike.append(float(option.log_strike)) price.append(float(option.price_in_forward_space)) ttm.append(float(option.ttm)) vol.append(float(option.implied_vol)) call_put.append(1 if option.option_type.is_call() else -1) return OptionArrays( options=options, - moneyness=np.array(moneyness), + log_strike=np.array(log_strike), price=np.array(price), ttm=np.array(ttm), implied_vol=np.array(vol), diff --git a/quantflow/sp/base.py b/quantflow/sp/base.py index e9c842ba..e60f3e94 100755 --- a/quantflow/sp/base.py +++ b/quantflow/sp/base.py @@ -103,6 +103,7 @@ class StochasticProcess1D(StochasticProcess): """ def marginal(self, t: FloatArrayLike) -> StochasticProcess1DMarginal: + """Marginal distribution of the process at time `t`""" return StochasticProcess1DMarginal(process=self, t=t) def domain_range(self) -> Bounds: @@ -167,8 +168,14 @@ def support(self, points: int = 100, *, std_mult: float = 4) -> FloatArray: float(self.mean()), std_mult * float(self.std()), points ) - def option_alpha(self) -> float: - """Option alpha parameter for integrability of call option transform""" + def call_option_carr_madan_alpha(self) -> float: + """Option alpha parameter for integrability of call option transform + in the Carr-Madan formula. + + 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. + """ return max(8 * np.max(np.exp(-2 * self.t)), 0.5) @@ -177,8 +184,12 @@ class IntensityProcess(StochasticProcess1D): as stochastic intensity """ - rate: float = Field(default=1.0, gt=0, description="Instantaneous initial rate") - kappa: float = Field(default=1.0, gt=0, description="Mean reversion speed") + rate: float = Field( + default=1.0, gt=0, description=r"Instantaneous initial rate $x_0$" + ) + kappa: float = Field( + default=1.0, gt=0, description=r"Mean reversion speed $\kappa$" + ) @abstractmethod def integrated_log_laplace( diff --git a/quantflow/sp/bns.py b/quantflow/sp/bns.py index 96139a29..a167e2bc 100644 --- a/quantflow/sp/bns.py +++ b/quantflow/sp/bns.py @@ -1,5 +1,7 @@ from __future__ import annotations +from typing import Self + import numpy as np from pydantic import Field from scipy.special import xlogy @@ -11,15 +13,42 @@ class BNS(StochasticProcess1D): - """Barndorff-Nielson--Shephard (BNS) stochastic volatility model""" + r"""Barndorff-Nielson & Shephard ([BNS](/bibliography#bns)) stochastic + volatility model. + + This is a stochastic volatility model where the variance process is given by a + non-Gaussian Ornstein-Uhlenbeck process driven by a pure-jump Lévy process. + The BNS model is defined by the following system of SDEs: + + \begin{equation} + \begin{aligned} + dx_t &= \sqrt{v_t} dw_t + \rho dz_{\kappa t} \\ + dv_t &= -\kappa v_t dt + dz_{\kappa t} + \end{aligned} + \end{equation} + + The model is flexible and can capture various stylized facts of financial markets, + such as volatility clustering and leverage effects. + + This implementation uses a [GammaOU][quantflow.sp.ou.GammaOU] process + for the variance, which is a common choice in the BNS model. + """ variance_process: GammaOU = Field( default_factory=GammaOU.create, description="Variance process" ) - rho: float = Field(default=0, ge=-1, le=1, description="Correlation") + rho: float = Field( + default=0, + gt=-1, + lt=1, + description="Correlation between variance and price processes, in (-1, 1)", + ) @classmethod - def create(cls, vol: float, kappa: float, decay: float, rho: float) -> BNS: + def create(cls, vol: float, kappa: float, decay: float, rho: float) -> Self: + """Convenience constructor for BNS process with parameters + of the variance process + """ return cls( variance_process=GammaOU.create(rate=vol * vol, kappa=kappa, decay=decay), rho=rho, diff --git a/quantflow/sp/ou.py b/quantflow/sp/ou.py index 61402168..908f0068 100755 --- a/quantflow/sp/ou.py +++ b/quantflow/sp/ou.py @@ -1,6 +1,7 @@ from __future__ import annotations -from typing import Generic +from abc import abstractmethod +from typing import Generic, Self import numpy as np from pydantic import Field @@ -23,21 +24,35 @@ class Vasicek(IntensityProcess): 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. - $$ - dx_t = \kappa (\theta - x_t) dt + \sigma dw_t - $$ + \begin{equation} + dx_t = \kappa (\theta - x_t) dt + \sigma dw_t + \end{equation} + + where $\kappa$ is the mean reversion speed, $\theta$ is the mean level, and + $\sigma$ is the volatility of the process. The solution to the SDE is given by + + \begin{equation} + x_t = x_0 e^{-\kappa t} + \theta (1 - e^{-\kappa t}) + + \sigma \int_0^t e^{-\kappa (t-s)} dw_s + \end{equation} """ bdlp: WeinerProcess = Field( default_factory=WeinerProcess, description="Background driving Weiner process", ) - theta: float = Field(default=1.0, gt=0, description="Mean rate") + 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 + + \begin{equation} + \phi_{x_t}(u) = - iu \mathbb{E}[x_t] + \frac{1}{2} u^2 \text{Var}[x_t] + \end{equation} + """ mu = self.analytical_mean(t) var = self.analytical_variance(t) - return u * (-complex(0, 1) * mu + 0.5 * var * u) + return u * (-1j * mu + 0.5 * var * u) def integrated_log_laplace(self, t: FloatArrayLike, u: Vector) -> Vector: raise NotImplementedError @@ -63,10 +78,28 @@ def domain_range(self) -> Bounds: return Bounds(-np.inf, np.inf) def analytical_mean(self, t: FloatArrayLike) -> FloatArrayLike: + r"""The analytical mean of the Vasicek process is given by + + \begin{equation} + \begin{aligned} + \mathbb{E}[x_t] &= x_0 e^{-\kappa t} + \theta (1 - e^{-\kappa t}) \\ + &= \theta \qquad \text{as } t \to \infty + \end{aligned} + \end{equation} + """ ekt = self.ekt(t) return self.rate * ekt + self.theta * (1 - ekt) def analytical_variance(self, t: FloatArrayLike) -> FloatArrayLike: + r"""The analytical variance of the Vasicek process is given by + + \begin{equation} + \begin{aligned} + \text{Var}[x_t] &= \frac{\sigma^2}{2\kappa} (1 - e^{-2\kappa t}) \\ + &= \frac{\sigma^2}{2\kappa} \qquad \text{as } t \to \infty + \end{aligned} + \end{equation} + """ ekt = self.ekt(2 * t) return 0.5 * self.bdlp.sigma2 * (1 - ekt) / self.kappa @@ -78,26 +111,66 @@ def analytical_cdf(self, t: FloatArrayLike, x: FloatArrayLike) -> FloatArrayLike class NGOU(IntensityProcess, Generic[D]): + r"""The general definition of a non-Gaussian OU process is defined as + the solution to the SDE + + \begin{equation} + dx_t = -\kappa x_t dt + dZ_{\kappa t} + \end{equation} + + where $Z_{\kappa t}$ is a pure jump Lévy process, also known as the background + driving Lévy process (BDLP). + The process is mean-reverting with mean reversion speed $\kappa>0$. + + The unusual timing $dZ_{\kappa t}$ is deliberately chosen so that it will turn + out that whatever the value of $\kappa$, the marginal distribution of $x_t$ + will be unchanged. + + The solution to the SDE is given by + + \begin{equation} + x_t = x_0 e^{-\kappa t} + \int_0^t e^{-\kappa (t-s)} dZ_{\kappa s} + \end{equation} + """ + bdlp: CompoundPoissonProcess[D] = Field( - description="Background driving Levy process", + description=( + "Background driving Lévy process is a " + "[CompoundPoissonProcess][quantflow.sp.poisson.CompoundPoissonProcess] " + "with jump distribution D" + ), ) + @property + def intensity(self) -> float: + """The intensity of the NGOU process is the intensity of the background + driving Lévy process""" + return self.bdlp.intensity + + @abstractmethod + def characteristic_exponent(self, t: FloatArrayLike, u: Vector) -> Vector: + r""" + \begin{equation} + \phi_{x_t,u} = iu x_0 e^{-\kappa t} + + \int_{ue^{-\kappa t}}^{u} \frac{\psi_Z(v)}{v} , dv + \end{equation} + """ + def sample_from_draws(self, draws: Paths, *args: Paths) -> Paths: raise NotImplementedError class GammaOU(NGOU[Exponential]): - @property - def intensity(self) -> float: - return self.bdlp.intensity + """Gamma OU process is a non-Gaussian OU process where the background + driving Lévy process is a compound Poisson process with exponential jumps""" @property def beta(self) -> float: - # TODO: find a better way for this return self.bdlp.jumps.decay @classmethod - def create(cls, rate: float = 1, decay: float = 1, kappa: float = 1) -> GammaOU: + def create(cls, rate: float = 1, decay: float = 1, kappa: float = 1) -> Self: + """Convenience constructor for GammaOU process""" return cls( rate=rate, kappa=kappa, @@ -107,6 +180,13 @@ def create(cls, rate: float = 1, decay: float = 1, kappa: float = 1) -> GammaOU: ) def characteristic_exponent(self, t: FloatArrayLike, u: Vector) -> Vector: + r""" + \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} + """ b = self.beta iu = Im * u c1 = iu * np.exp(-self.kappa * t) diff --git a/quantflow/utils/marginal.py b/quantflow/utils/marginal.py index 3a37dc6c..14b79d2c 100644 --- a/quantflow/utils/marginal.py +++ b/quantflow/utils/marginal.py @@ -1,11 +1,12 @@ from __future__ import annotations +import enum from abc import ABC, abstractmethod from typing import Callable, cast import numpy as np import pandas as pd -from pydantic import BaseModel +from pydantic import BaseModel, Field from scipy.optimize import Bounds from typing_extensions import Annotated, Doc @@ -13,70 +14,329 @@ from .types import FloatArray, FloatArrayLike, Vector -class Marginal1D(BaseModel, ABC, extra="forbid"): - """Marginal distribution""" +class OptionPricingMethod(enum.StrEnum): + """Method to use for option pricing via Fourier transform of the call option + price""" - @abstractmethod - def characteristic(self, u: Vector) -> Vector: - """ - Compute the characteristic function on support points `n`. - """ + CARR_MADAN = enum.auto() + LEWIS = enum.auto() + COS = enum.auto() + + +class OptionPricingResult(BaseModel, ABC): + method: OptionPricingMethod @abstractmethod - def support(self, points: int = 100, *, std_mult: float = 3) -> FloatArray: - """ - Compute the x axis. - """ + def call_at(self, log_strikes: FloatArrayLike) -> FloatArray: + """Evaluate call prices at arbitrary log-strikes.""" - def mean(self) -> FloatArrayLike: - """Expected value - This should be overloaded if a more efficient way of computing the mean - """ - return self.mean_from_characteristic() +class OptionPricingTransformResult(OptionPricingResult, arbitrary_types_allowed=True): + log_strikes: FloatArray + call_prices: FloatArray - def variance(self) -> FloatArrayLike: - """Variance + def call_at(self, log_strikes: FloatArrayLike) -> FloatArray: + return cast( + FloatArray, + np.interp(log_strikes, self.log_strikes, self.call_prices), + ) - This should be overloaded if a more efficient way of computing the - """ - return self.variance_from_characteristic() - def domain_range(self) -> Bounds: - """The space domain range for the random variable +class OptionPricingCosResult(OptionPricingResult, arbitrary_types_allowed=True): + """Result of call option pricing via the COS method. - This should be overloaded if required - """ - return default_bounds() + Stores the precomputed coefficient vector, enabling $O(N)$ evaluation + at any log-strike without recomputing the characteristic function. + """ - def frequency_range(self, max_frequency: float | None = None) -> float: - """The frequency domain range for the characteristic function + a: float = Field(description="Left endpoint of the truncation interval [a, b]") + nu: FloatArray = Field(description="Cosine frequency grid j*pi/(b-a)") + coeff: np.ndarray = Field( + description="Complex coefficient vector w_j * phi(nu_j) * V_j" + ) - This should be overloaded if required + def call_at( + self, + log_strikes: Annotated[ + FloatArrayLike, Doc("Log-strikes at which to evaluate the call price") + ], + ) -> FloatArray: + """Evaluate call prices at arbitrary log-strikes in $O(N)$ per strike.""" + k = np.asarray(log_strikes) + return np.real( + np.exp(-1j * np.outer(k + self.a, self.nu)) @ self.coeff + ) * np.exp(k) + + +class Marginal1D(BaseModel, ABC, extra="forbid"): + r"""Abstract 1D marginal distribution with Fourier-based option pricing. + + This class represents the marginal distribution. + It provides methods to compute the pdf, cdf, and option + prices from the characteristic function of the underlying process, as well as + their Jacobians with respect to the parameters of the process. + + Option prices are computed via Fourier inversion of the + [characteristic function](/glossary#characteristic-function), + using two supported formulas: + + * [Carr & Madan](/bibliography#carr_madan): uses a damping parameter $\alpha$ + to ensure integrability. + The integrand is evaluated along a contour shifted by $\alpha$ in the + imaginary direction. See [call_option_carr_madan][.call_option_carr_madan]. + * [Lewis](/bibliography#lewis): no damping parameter required. + The contour is fixed at + imaginary part $1/2$, giving an integrand that is naturally bounded for + all real $u$. See [call_option_lewis][.call_option_lewis]. + * [COS method](/bibliography#cos): uses a Fourier-cosine expansion of the + option payoff, with coefficients that depend on the characteristic function + evaluated at a cosine frequency grid. See [call_option_cos][.call_option_cos]. + + It is the base class for the + [StochasticProcess1DMarginal][quantflow.sp.base.StochasticProcess1DMarginal]. + """ + + def call_option( + self, + n: Annotated[ + int | None, + Doc("Number of discretization points for the transform. Defaults to 128."), + ] = None, + *, + pricing_method: Annotated[ + OptionPricingMethod, + Doc( + "Method to use for option pricing via Fourier transform of the " + "call option price. Defaults to Lewis." + ), + ] = OptionPricingMethod.CARR_MADAN, + max_frequency: Annotated[ + float | None, + Doc( + "Maximum frequency for the transform grid. " + "Defaults to frequency_range()." + ), + ] = None, + max_log_strike: Annotated[ + float, Doc("Maximum absolute log-strike for the output grid.") + ] = 1, + alpha: Annotated[ + float | None, + Doc( + "Damping parameter for integrability of the Carr-Madan integrand, " + "it is ignored if not applicable." + ), + ] = None, + simpson_rule: Annotated[ + bool, Doc("Use Simpson's rule for integration. Default is False.") + ] = False, + use_fft: Annotated[ + bool, Doc("Use FFT for the transform rather than FRFT. Default is False.") + ] = False, + ) -> OptionPricingResult: + """Price call options via one of the available Fourier-based methods, + as a function of log-strike """ - return Bounds(0, max_frequency or 20) + match pricing_method: + case OptionPricingMethod.COS: + return self.call_option_cos( + n, + max_log_strike=max_log_strike, + ) + case OptionPricingMethod.CARR_MADAN: + return self.call_option_carr_madan( + n, + max_frequency=max_frequency, + max_log_strike=max_log_strike, + alpha=alpha, + simpson_rule=simpson_rule, + use_fft=use_fft, + ) + case _: + return self.call_option_lewis( + n, + max_frequency=max_frequency, + max_log_strike=max_log_strike, + simpson_rule=simpson_rule, + use_fft=use_fft, + ) + + def call_option_carr_madan( + self, + n: Annotated[ + int | None, + Doc("Number of discretization points for the transform. Defaults to 128."), + ] = None, + *, + max_frequency: Annotated[ + float | None, + Doc( + "Maximum frequency for the transform grid. " + "Defaults to frequency_range()." + ), + ] = None, + max_log_strike: Annotated[ + float, Doc("Maximum absolute log-strike for the output grid.") + ] = 1, + alpha: Annotated[ + float | None, + Doc( + "Damping parameter for integrability of the Carr-Madan integrand. " + "Defaults to " + "[call_option_carr_madan_alpha][..call_option_carr_madan_alpha]." + ), + ] = None, + simpson_rule: Annotated[ + bool, Doc("Use Simpson's rule for integration. Default is False.") + ] = False, + use_fft: Annotated[ + bool, Doc("Use FFT for the transform. Default is False.") + ] = False, + ) -> OptionPricingTransformResult: + """Call option price via Carr & Madan method""" + transform = self.get_transform( + n, + lambda m: self.option_support(m + 1, max_log_strike=max_log_strike), + max_frequency=max_frequency, + simpson_rule=simpson_rule, + use_fft=use_fft, + ) + alpha = alpha or self.call_option_carr_madan_alpha() + phi = cast( + np.ndarray, + self.call_option_transform(transform.frequency_domain - 1j * alpha), + ) + result = transform(phi, use_fft=use_fft) + return OptionPricingTransformResult( + method=OptionPricingMethod.CARR_MADAN, + log_strikes=result.x, + call_prices=result.y * np.exp(-alpha * result.x), + ) - def std(self) -> FloatArrayLike: - """Standard deviation at a time horizon""" - return np.sqrt(self.variance()) + 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 - def mean_from_characteristic(self, *, d: float = 0.001) -> FloatArrayLike: - """Calculate mean as first derivative of characteristic function at 0""" - m = -0.5 * 1j * (self.characteristic(d) - self.characteristic(-d)) / d - return m.real + def call_option_cos( + self, + n: Annotated[ + int | None, + Doc("Number of cosine series terms. Defaults to 128."), + ] = None, + *, + max_log_strike: Annotated[ + float, Doc("Maximum absolute log-strike for the output grid.") + ] = 1, + L: Annotated[ + float, + Doc( + "Truncation parameter: the integration interval is set to " + "[-L*std, L*std]." + ), + ] = 12.0, + ) -> OptionPricingCosResult: + r"""Call option price via the COS method (Fang & Oosterlee 2008). + + The call price at log-strike $k$ is approximated by the cosine series + + \begin{equation} + C(k) \approx e^k \sum_{j=0}^{N-1}{}^{\prime} + \mathrm{Re}\!\left[ + \phi\!\left(\frac{j\pi}{b-a}\right) + e^{-i j\pi (k+a)/(b-a)} + \right] V_j + \end{equation} + + where the prime denotes a half weight on the $j=0$ term, $[a,b]$ is the + truncation interval, and $V_j$ are the cosine payoff coefficients for the + normalised call payoff $(e^y - 1)^+$ integrated over $[0, b]$. + The $e^k$ factor converts from strike-normalised to forward-space pricing. + + Returns an [OptionPricingCosResult][.OptionPricingCosResult] with the + precomputed coefficient vector. Use + [call_at][quantflow.utils.marginal.OptionPricingCosResult.call_at] + to evaluate at arbitrary log-strikes in $O(N)$ per strike. + """ + n = n or 128 + std = float(self.std()) + a = -L * std + b = L * std + bma = b - a + + j = np.arange(n) + nu = j * np.pi / bma + + phi = np.asarray(self.characteristic_corrected(nu)) + + sin_nu_a = np.sin(nu * a) + cos_nu_a = np.cos(nu * a) + sign_j = (-1.0) ** j + safe_nu = np.where(j == 0, 1.0, nu) + chi = (np.exp(b) * sign_j - cos_nu_a + nu * sin_nu_a) / np.where( + j == 0, 1.0, 1.0 + nu**2 + ) + psi = np.where(j == 0, b, sin_nu_a / safe_nu) + V = 2.0 / bma * (chi - psi) + weights = np.ones(n) + weights[0] = 0.5 + return OptionPricingCosResult( + method=OptionPricingMethod.COS, + a=a, + nu=nu, + coeff=weights * phi * V, + ) - def std_from_characteristic(self) -> FloatArrayLike: - """Calculate standard deviation as square root of variance""" - return np.sqrt(self.variance_from_characteristic()) + def call_option_lewis( + self, + n: Annotated[ + int | None, + Doc("Number of discretization points for the transform. Defaults to 128."), + ] = None, + *, + max_frequency: Annotated[ + float | None, + Doc( + "Maximum frequency for the transform grid. " + "Defaults to frequency_range()." + ), + ] = None, + max_log_strike: Annotated[ + float, Doc("Maximum absolute log-strike for the output grid.") + ] = 1, + simpson_rule: Annotated[ + bool, Doc("Use Simpson's rule for integration. Default is False.") + ] = False, + use_fft: Annotated[ + bool, Doc("Use FFT for the transform. Default is False.") + ] = False, + ) -> OptionPricingTransformResult: + """Call option price via the Lewis (2001) formula""" + transform = self.get_transform( + n, + lambda m: self.option_support(m + 1, max_log_strike=max_log_strike), + max_frequency=max_frequency, + simpson_rule=simpson_rule, + use_fft=use_fft, + ) + phi = cast(np.ndarray, self.lewis_transform(transform.frequency_domain)) + result = transform(phi, use_fft=use_fft) + k = result.x + return OptionPricingTransformResult( + log_strikes=k, + call_prices=1.0 - np.exp(0.5 * k) * result.y, + method=OptionPricingMethod.LEWIS, + ) - def variance_from_characteristic(self, *, d: float = 0.001) -> FloatArrayLike: - """Calculate variance as second derivative of characteristic function at 0""" - c1 = self.characteristic(d) - c0 = self.characteristic(0) - c2 = self.characteristic(-d) - m = -0.5 * 1j * (c1 - c2) / d - s = -(c1 - 2 * c0 + c2) / (d * d) - m * m - return s.real + def call_option_transform( + self, + u: Annotated[ + Vector, Doc("Frequency domain points (possibly complex-shifted).") + ], + ) -> Vector: + """Call option transform""" + uj = 1j * u + return self.characteristic_corrected(u - 1j) / (uj * uj + uj) def cdf( self, @@ -93,116 +353,225 @@ def cdf( """ raise NotImplementedError("Analytical CDF not available") - def pdf( + def cdf_from_characteristic( self, - x: Annotated[ - FloatArrayLike, + n: Annotated[ + int | None, + Doc("Number of discretization points for the transform. Defaults to 128."), + ] = None, + *, + max_frequency: Annotated[ + float | None, Doc( - "Location in the stochastic process domain space. If a numpy array," - " the output should have the same shape as the input." + "Maximum frequency for the transform grid. " + "Defaults to frequency_range()." ), + ] = None, + simpson_rule: Annotated[ + bool, Doc("Use Simpson's rule for integration. Default is False.") + ] = False, + use_fft: Annotated[ + bool, Doc("Use FFT for the transform. Default is False.") + ] = False, + frequency_n: Annotated[ + int | None, + Doc("Number of points for the frequency grid. Overrides n if provided."), + ] = None, + ) -> TransformResult: + raise NotImplementedError("CDF not available") + + def cdf_jacobian( + self, + x: Annotated[ + FloatArrayLike, Doc("Location in the state space of the process.") ], - ) -> FloatArrayLike: + ) -> np.ndarray: """ - Computes the probability density (or mass) function of the process. + Jacobian of the cdf with respect to the parameters of the process. + It is useful for optimization purposes if necessary. - It has a base implementation that computes the pdf from the - [cdf][quantflow.utils.marginal.Marginal1D.cdf] method, but a subclass should - overload this method if a - more optimized way of computing it is available. + Optional to implement, otherwise raises ``NotImplementedError`` if called. """ - raise NotImplementedError("Analytical PDF not available") + raise NotImplementedError("Analytical CDF Jacobian not available") - def pdf_from_characteristic( + @abstractmethod + def characteristic( + self, + u: Annotated[Vector, Doc("Frequency domain points.")], + ) -> Vector: + """ + Compute the characteristic function on frequency domain points $u$ + """ + + def characteristic_corrected( + self, + u: Annotated[Vector, Doc("Frequency domain points.")], + ) -> Vector: + """Characteristic function corrected for the convexity of the log-price + distribution""" + convexity = np.log(self.characteristic(-1j)) + return self.characteristic(u) * np.exp(-1j * u * convexity) + + def characteristic_df( self, n: Annotated[ int | None, - Doc( - "Number of discretization points to use in the transform." - " If None, use 128." - ), + Doc("Number of discretization points for the transform. Defaults to 128."), ] = None, *, max_frequency: Annotated[ float | None, Doc( - "The maximum frequency to use in the transform. If not provided," - " the value from the [frequency_range]" - "[quantflow.utils.marginal.Marginal1D.frequency_range] method is used." - " Only needed for special cases/testing." + "Maximum frequency for the transform grid. " + "Defaults to frequency_range()." ), ] = None, simpson_rule: Annotated[ bool, Doc("Use Simpson's rule for integration. Default is False.") ] = False, - use_fft: Annotated[ - bool, Doc("Use FFT for the transform. Default is False.") - ] = False, - frequency_n: int | None = None, - ) -> TransformResult: + ) -> pd.DataFrame: """ - Compute the probability density function from the characteristic function. + Compute the characteristic function with n discretization points + and a max frequency """ transform = self.get_transform( - frequency_n or n, + n, self.support, max_frequency=max_frequency, simpson_rule=simpson_rule, - use_fft=use_fft, ) - psi = cast(np.ndarray, self.characteristic(transform.frequency_domain)) - return transform(psi, use_fft=use_fft) + psi = self.characteristic(transform.frequency_domain) + return transform.characteristic_df(cast(np.ndarray, psi)) - def cdf_from_characteristic( + def domain_range(self) -> Bounds: + """The space domain range for the random variable + + This should be overloaded if required + """ + return default_bounds() + + def frequency_range( self, - n: int | None = None, - *, - max_frequency: float | None = None, - simpson_rule: bool = False, - use_fft: bool = False, - frequency_n: int | None = None, - ) -> TransformResult: - raise NotImplementedError("CDF not available") + max_frequency: Annotated[ + float | None, + Doc("Upper bound of the frequency grid. Defaults to 20 if None."), + ] = None, + ) -> float: + """The frequency domain range for the characteristic function - def call_option( + This should be overloaded if required + """ + return Bounds(0, max_frequency or 20) + + def get_transform( self, - n: int | None = None, + n: Annotated[ + int | None, + Doc("Number of discretization points. Defaults to 128."), + ], + support: Annotated[ + Callable[[int], FloatArray], + Doc("Function returning the space domain grid given the number of points."), + ], *, - max_frequency: float | None = None, - max_moneyness: float = 1, - alpha: float | None = None, - simpson_rule: bool = False, - use_fft: bool = False, - ) -> TransformResult: - transform = self.get_transform( + max_frequency: Annotated[ + float | None, + Doc( + "Maximum frequency for the transform grid. " + "Defaults to frequency_range()." + ), + ] = None, + simpson_rule: Annotated[ + bool, Doc("Use Simpson's rule for integration. Default is False.") + ] = False, + use_fft: Annotated[ + bool, Doc("Use FFT for the transform. Default is False.") + ] = False, + ) -> Transform: + n = n or 128 + if use_fft: + bounds = self.domain_range() + else: + x = support(n) + bounds = Bounds(float(np.min(x)), float(np.max(x))) + return Transform.create( n, - lambda m: self.option_support(m + 1, max_moneyness=max_moneyness), - max_frequency=max_frequency, + frequency_range=self.frequency_range(max_frequency), + domain_range=bounds, simpson_rule=simpson_rule, - use_fft=use_fft, - ) - alpha = alpha or self.option_alpha() - phi = cast( - np.ndarray, - self.call_option_transform(transform.frequency_domain - 1j * alpha), ) - result = transform(phi, use_fft=use_fft) - return TransformResult(x=result.x, y=result.y * np.exp(-alpha * result.x)) + + def lewis_transform( + self, + u: Annotated[Vector, Doc("Frequency domain points.")], + ) -> Vector: + """Lewis (2001) call option transform - no damping parameter required""" + return self.characteristic_corrected(u - 0.5j) / (u * u + 0.25) + + def mean(self) -> FloatArrayLike: + """Expected value + + By default it uses the + [mean_from_characteristic][..mean_from_characteristic] method. + This should be overloaded if a more efficient/analytical way of computing + the mean is available. + """ + return self.mean_from_characteristic() + + def mean_from_characteristic( + self, + *, + d: Annotated[ + float, + Doc("Step size for finite-difference approximation of the derivative."), + ] = 0.001, + ) -> FloatArrayLike: + """Calculate mean as first derivative of characteristic function at 0""" + m = -0.5 * 1j * (self.characteristic(d) - self.characteristic(-d)) / d + return m.real + + def option_support( + self, + points: Annotated[int, Doc("Number of support points.")] = 101, + max_log_strike: Annotated[float, Doc("Maximum absolute log-strike.")] = 1.0, + ) -> FloatArray: + """ + Compute the x axis. + """ + return np.linspace(-max_log_strike, max_log_strike, points) def option_time_value( self, - n: int = 128, + n: Annotated[ + int, + Doc("Number of discretization points for the transform."), + ] = 128, *, - max_frequency: float | None = None, - max_moneyness: float = 1, - alpha: float = 1.1, - simpson_rule: bool = False, - use_fft: bool = False, + max_frequency: Annotated[ + float | None, + Doc( + "Maximum frequency for the transform grid. " + "Defaults to frequency_range()." + ), + ] = None, + max_log_strike: Annotated[ + float, Doc("Maximum absolute log-strike for the output grid.") + ] = 1, + alpha: Annotated[ + float, + Doc("Contour shift parameter controlling the integration strip."), + ] = 1.1, + simpson_rule: Annotated[ + bool, Doc("Use Simpson's rule for integration. Default is False.") + ] = False, + use_fft: Annotated[ + bool, Doc("Use FFT for the transform. Default is False.") + ] = False, ) -> TransformResult: """Option time value""" transform = self.get_transform( n, - lambda m: self.option_support(m + 1, max_moneyness=max_moneyness), + lambda m: self.option_support(m + 1, max_log_strike=max_log_strike), max_frequency=max_frequency, simpson_rule=simpson_rule, use_fft=use_fft, @@ -215,49 +584,92 @@ def option_time_value( time_value = result.y / np.sinh(alpha * result.x) return TransformResult(x=result.x, y=time_value) - def characteristic_df( + def option_time_value_transform( + self, + u: Annotated[Vector, Doc("Frequency domain points.")], + alpha: Annotated[ + float, Doc("Contour shift parameter controlling the integration strip.") + ] = 1.1, + ) -> Vector: + """Option time value transform + + This transform does not require any additional correction since + the integrand is already bounded for positive and negative moneyness""" + ia = 1j * alpha + return 0.5 * ( + self._option_time_value_transform(u - ia) + - self._option_time_value_transform(u + ia) + ) + + def pdf( + self, + x: Annotated[ + FloatArrayLike, + Doc( + "Location in the stochastic process domain space. If a numpy array," + " the output should have the same shape as the input." + ), + ], + ) -> FloatArrayLike: + """ + Computes the probability density (or mass) function of the process. + + It has a base implementation that computes the pdf from the + [cdf][quantflow.utils.marginal.Marginal1D.cdf] method, but a subclass should + overload this method if a + more optimized way of computing it is available. + """ + raise NotImplementedError("Analytical PDF not available") + + def pdf_from_characteristic( self, - n: int | None = None, + n: Annotated[ + int | None, + Doc( + "Number of discretization points to use in the transform." + " If None, use 128." + ), + ] = None, *, - max_frequency: float | None = None, - simpson_rule: bool = False, - ) -> pd.DataFrame: + max_frequency: Annotated[ + float | None, + Doc( + "The maximum frequency to use in the transform. If not provided," + " the value from the [frequency_range]" + "[quantflow.utils.marginal.Marginal1D.frequency_range] method is used." + " Only needed for special cases/testing." + ), + ] = None, + simpson_rule: Annotated[ + bool, Doc("Use Simpson's rule for integration. Default is False.") + ] = False, + use_fft: Annotated[ + bool, Doc("Use FFT for the transform. Default is False.") + ] = False, + frequency_n: Annotated[ + int | None, + Doc("Number of points for the frequency grid. Overrides n if provided."), + ] = None, + ) -> TransformResult: """ - Compute the characteristic function with n discretization points - and a max frequency + Compute the probability density function from the characteristic function. """ transform = self.get_transform( - n, + frequency_n or n, self.support, max_frequency=max_frequency, simpson_rule=simpson_rule, + use_fft=use_fft, ) - psi = self.characteristic(transform.frequency_domain) - return transform.characteristic_df(cast(np.ndarray, psi)) + psi = cast(np.ndarray, self.characteristic(transform.frequency_domain)) + return transform(psi, use_fft=use_fft) - def get_transform( + def pdf_jacobian( self, - n: int | None, - support: Callable[[int], FloatArray], - *, - max_frequency: float | None = None, - simpson_rule: bool = False, - use_fft: bool = False, - ) -> Transform: - n = n or 128 - if use_fft: - bounds = self.domain_range() - else: - x = support(n) - bounds = Bounds(float(np.min(x)), float(np.max(x))) - return Transform.create( - n, - frequency_range=self.frequency_range(max_frequency), - domain_range=bounds, - simpson_rule=simpson_rule, - ) - - def pdf_jacobian(self, x: FloatArrayLike) -> FloatArrayLike: + x: Annotated[ + FloatArrayLike, Doc("Location in the state space of the process.") + ], + ) -> FloatArrayLike: """ Jacobian of the pdf with respect to the parameters of the process. It has a base implementation that computes it from the @@ -267,46 +679,57 @@ def pdf_jacobian(self, x: FloatArrayLike) -> FloatArrayLike: """ return self.cdf_jacobian(x) - self.cdf_jacobian(x - 1) - def cdf_jacobian(self, x: FloatArrayLike) -> np.ndarray: - """ - Jacobian of the cdf with respect to the parameters of the process. - It is useful for optimization purposes if necessary. + def std(self) -> FloatArrayLike: + """Standard deviation at a time horizon""" + return np.sqrt(self.variance()) - Optional to implement, otherwise raises ``NotImplementedError`` if called. - """ - raise NotImplementedError("Analytical CDF Jacobian not available") + def std_from_characteristic(self) -> FloatArrayLike: + """Calculate standard deviation as square root of variance""" + return np.sqrt(self.variance_from_characteristic()) - def option_support( - self, points: int = 101, max_moneyness: float = 1.0 + @abstractmethod + def support( + self, + points: Annotated[int, Doc("Number of support points.")] = 100, + *, + std_mult: Annotated[ + float, Doc("Standard deviation multiplier for the support range.") + ] = 3, ) -> FloatArray: """ Compute the x axis. """ - return np.linspace(-max_moneyness, max_moneyness, points) - - # Fourier Transforms for options - - def call_option_transform(self, u: Vector) -> Vector: - """Call option transform""" - uj = 1j * u - return self.characteristic_corrected(u - 1j) / (uj * uj + uj) - def characteristic_corrected(self, u: Vector) -> Vector: - convexity = np.log(self.characteristic(-1j)) - return self.characteristic(u) * np.exp(-1j * u * convexity) + def variance(self) -> FloatArrayLike: + """Variance - def option_time_value_transform(self, u: Vector, alpha: float = 1.1) -> Vector: - """Option time value transform + By default it uses the + [variance_from_characteristic][..variance_from_characteristic] method. + This should be overloaded if a more efficient/analytical way of computing + the variance is available. + """ + return self.variance_from_characteristic() - This transform does not require any additional correction since - the integrand is already bounded for positive and negative moneyness""" - ia = 1j * alpha - return 0.5 * ( - self._option_time_value_transform(u - ia) - - self._option_time_value_transform(u + ia) - ) + def variance_from_characteristic( + self, + *, + d: Annotated[ + float, + Doc("Step size for finite-difference approximation of the derivative."), + ] = 0.001, + ) -> FloatArrayLike: + """Calculate variance as second derivative of characteristic function at 0""" + c1 = self.characteristic(d) + c0 = self.characteristic(0) + c2 = self.characteristic(-d) + m = -0.5 * 1j * (c1 - c2) / d + s = -(c1 - 2 * c0 + c2) / (d * d) - m * m + return s.real - def _option_time_value_transform(self, u: Vector) -> Vector: + def _option_time_value_transform( + self, + u: Annotated[Vector, Doc("Frequency domain points (complex-shifted).")], + ) -> Vector: """Option time value transform This transform does not require any additional correction since @@ -315,7 +738,3 @@ def _option_time_value_transform(self, u: Vector) -> Vector: return ( 1 / (1 + iu) - 1 / iu - self.characteristic_corrected(u - 1j) / (u * u - iu) ) - - def option_alpha(self) -> float: - """Option alpha""" - return 2.0 diff --git a/quantflow/utils/plot.py b/quantflow/utils/plot.py index b11dfc77..8e8eb88c 100644 --- a/quantflow/utils/plot.py +++ b/quantflow/utils/plot.py @@ -106,7 +106,7 @@ def plot_vol_surface( *, model: pd.DataFrame | None = None, marker_size: int = 10, - x_series: str = "moneyness_ttm", + x_series: str = "moneyness", series: str = "implied_vol", color_series: str = "side", fig: Any | None = None, @@ -135,7 +135,7 @@ def plot_vol_surface( if model is not None: fig_.add_trace( go.Scatter( - x=model["moneyness_ttm"], + x=model["moneyness"], y=model[series], name="model", mode="lines", @@ -156,7 +156,7 @@ def plot_vol_surface_3d( **kwargs: Any, ) -> Any: check_plotly() - new_fig = px.scatter_3d(df, x="moneyness_ttm", y="ttm", z=series, color="side") + new_fig = px.scatter_3d(df, x="moneyness", y="ttm", z=series, color="side") if fig is None: fig = new_fig fig.update_layout(scene_dragmode=dragmode, uirevision=uirevision, **kwargs) @@ -179,7 +179,7 @@ def plot_vol_cross( fig = fig or go.Figure() fig.add_trace( go.Scatter( - x=data["moneyness_ttm"], + x=data["moneyness"], y=data[series], name=name, mode="lines", @@ -188,13 +188,13 @@ def plot_vol_cross( if data2 is not None: fig.add_trace( go.Scatter( - x=data2["moneyness_ttm"], + x=data2["moneyness"], y=data2[series], name="model", mode="lines", ) ) - return fig.update_layout(xaxis_title="moneyness_ttm", yaxis_title=series) + return fig.update_layout(xaxis_title="moneyness", yaxis_title=series) def plot3d( diff --git a/quantflow_tests/test_options.py b/quantflow_tests/test_options.py index bf7dc225..c35eb9ad 100644 --- a/quantflow_tests/test_options.py +++ b/quantflow_tests/test_options.py @@ -168,7 +168,7 @@ def test_black_vol(vol_surface: VolSurface): def test_call_put_parity(): option = OptionPrice.create(100).calculate_price() - assert option.moneyness == 0 + assert option.log_strike == 0 assert option.price == option.call_price option2 = OptionPrice.create(100, option_type=OptionType.put).calculate_price() assert option2.price == option2.put_price @@ -178,7 +178,7 @@ def test_call_put_parity(): def test_call_put_parity_otm(): option = OptionPrice.create(105, forward=100).calculate_price() - assert option.moneyness > 0 + assert option.log_strike > 0 assert option.price == option.call_price option2 = OptionPrice.create( 105, forward=100, option_type=OptionType.put @@ -231,7 +231,7 @@ def _synthetic_options( ref_date=ref_date, maturity=maturity, ) - entry = OptionEntry(ttm=ttm, moneyness=float(k)) + entry = OptionEntry(ttm=ttm, log_strike=float(k)) entry.options.append(op) options[key] = entry return options diff --git a/uv.lock b/uv.lock index e7914b16..f654d197 100644 --- a/uv.lock +++ b/uv.lock @@ -712,7 +712,7 @@ name = "cuda-bindings" version = "12.9.6" source = { registry = "https://pypi.org/simple" } dependencies = [ - { name = "cuda-pathfinder" }, + { name = "cuda-pathfinder", marker = "sys_platform != 'emscripten' and sys_platform != 'win32'" }, ] wheels = [ { url = "https://files.pythonhosted.org/packages/1f/a5/e9d37c10f6c27c9c65d53c6cd6d9763e1df99c004780585fc2ad9041fbe3/cuda_bindings-12.9.6-cp311-cp311-manylinux_2_24_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:2662f59db67d9aeaf8959c593c91f600792c2970cf02cae2814387fc687b115a", size = 7090971, upload-time = "2026-03-11T14:47:29.526Z" }, @@ -747,37 +747,37 @@ wheels = [ [package.optional-dependencies] cublas = [ - { name = "nvidia-cublas-cu12", marker = "sys_platform == 'linux' or sys_platform == 'win32'" }, + { name = "nvidia-cublas-cu12", marker = "sys_platform == 'linux'" }, ] cudart = [ - { name = "nvidia-cuda-runtime-cu12", marker = "sys_platform == 'linux' or sys_platform == 'win32'" }, + { name = "nvidia-cuda-runtime-cu12", marker = "sys_platform == 'linux'" }, ] cufft = [ - { name = "nvidia-cufft-cu12", marker = "sys_platform == 'linux' or sys_platform == 'win32'" }, + { name = "nvidia-cufft-cu12", marker = "sys_platform == 'linux'" }, ] cufile = [ { name = "nvidia-cufile-cu12", marker = "sys_platform == 'linux'" }, ] cupti = [ - { name = "nvidia-cuda-cupti-cu12", marker = "sys_platform == 'linux' or sys_platform == 'win32'" }, + { name = "nvidia-cuda-cupti-cu12", marker = "sys_platform == 'linux'" }, ] curand = [ - { name = "nvidia-curand-cu12", marker = "sys_platform == 'linux' or sys_platform == 'win32'" }, + { name = "nvidia-curand-cu12", marker = "sys_platform == 'linux'" }, ] cusolver = [ - { name = "nvidia-cusolver-cu12", marker = "sys_platform == 'linux' or sys_platform == 'win32'" }, + { name = "nvidia-cusolver-cu12", marker = "sys_platform == 'linux'" }, ] cusparse = [ - { name = "nvidia-cusparse-cu12", marker = "sys_platform == 'linux' or sys_platform == 'win32'" }, + { name = "nvidia-cusparse-cu12", marker = "sys_platform == 'linux'" }, ] nvjitlink = [ - { name = "nvidia-nvjitlink-cu12", marker = "sys_platform == 'linux' or sys_platform == 'win32'" }, + { name = "nvidia-nvjitlink-cu12", marker = "sys_platform == 'linux'" }, ] nvrtc = [ - { name = "nvidia-cuda-nvrtc-cu12", marker = "sys_platform == 'linux' or sys_platform == 'win32'" }, + { name = "nvidia-cuda-nvrtc-cu12", marker = "sys_platform == 'linux'" }, ] nvtx = [ - { name = "nvidia-nvtx-cu12", marker = "sys_platform == 'linux' or sys_platform == 'win32'" }, + { name = "nvidia-nvtx-cu12", marker = "sys_platform == 'linux'" }, ] [[package]] @@ -2322,7 +2322,7 @@ name = "nvidia-cudnn-cu12" version = "9.10.2.21" source = { registry = "https://pypi.org/simple" } dependencies = [ - { name = "nvidia-cublas-cu12" }, + { name = "nvidia-cublas-cu12", marker = "sys_platform != 'emscripten' and sys_platform != 'win32'" }, ] wheels = [ { url = "https://files.pythonhosted.org/packages/fa/41/e79269ce215c857c935fd86bcfe91a451a584dfc27f1e068f568b9ad1ab7/nvidia_cudnn_cu12-9.10.2.21-py3-none-manylinux_2_27_aarch64.whl", hash = "sha256:c9132cc3f8958447b4910a1720036d9eff5928cc3179b0a51fb6d167c6cc87d8", size = 705026878, upload-time = "2025-06-06T21:52:51.348Z" }, @@ -2334,7 +2334,7 @@ name = "nvidia-cufft-cu12" version = "11.3.0.4" source = { registry = "https://pypi.org/simple" } dependencies = [ - { name = "nvidia-nvjitlink-cu12" }, + { name = "nvidia-nvjitlink-cu12", marker = "sys_platform != 'emscripten' and sys_platform != 'win32'" }, ] wheels = [ { url = "https://files.pythonhosted.org/packages/1f/37/c50d2b2f2c07e146776389e3080f4faf70bcc4fa6e19d65bb54ca174ebc3/nvidia_cufft_cu12-11.3.0.4-py3-none-manylinux2014_aarch64.manylinux_2_17_aarch64.whl", hash = "sha256:d16079550df460376455cba121db6564089176d9bac9e4f360493ca4741b22a6", size = 200164144, upload-time = "2024-11-20T17:40:58.288Z" }, @@ -2368,9 +2368,9 @@ name = "nvidia-cusolver-cu12" version = "11.7.1.2" source = { registry = "https://pypi.org/simple" } dependencies = [ - { name = "nvidia-cublas-cu12" }, - { name = "nvidia-cusparse-cu12" }, - { name = "nvidia-nvjitlink-cu12" }, + { name = "nvidia-cublas-cu12", marker = "sys_platform != 'emscripten' and sys_platform != 'win32'" }, + { name = "nvidia-cusparse-cu12", marker = "sys_platform != 'emscripten' and sys_platform != 'win32'" }, + { name = "nvidia-nvjitlink-cu12", marker = "sys_platform != 'emscripten' and sys_platform != 'win32'" }, ] wheels = [ { url = "https://files.pythonhosted.org/packages/93/17/dbe1aa865e4fdc7b6d4d0dd308fdd5aaab60f939abfc0ea1954eac4fb113/nvidia_cusolver_cu12-11.7.1.2-py3-none-manylinux2014_aarch64.whl", hash = "sha256:0ce237ef60acde1efc457335a2ddadfd7610b892d94efee7b776c64bb1cac9e0", size = 157833628, upload-time = "2024-10-01T17:05:05.591Z" }, @@ -2384,7 +2384,7 @@ name = "nvidia-cusparse-cu12" version = "12.5.4.2" source = { registry = "https://pypi.org/simple" } dependencies = [ - { name = "nvidia-nvjitlink-cu12" }, + { name = "nvidia-nvjitlink-cu12", marker = "sys_platform != 'emscripten' and sys_platform != 'win32'" }, ] wheels = [ { url = "https://files.pythonhosted.org/packages/eb/eb/6681efd0aa7df96b4f8067b3ce7246833dd36830bb4cec8896182773db7d/nvidia_cusparse_cu12-12.5.4.2-py3-none-manylinux2014_aarch64.manylinux_2_17_aarch64.whl", hash = "sha256:d25b62fb18751758fe3c93a4a08eff08effedfe4edf1c6bb5afd0890fe88f887", size = 216451147, upload-time = "2024-11-20T17:44:18.055Z" }, @@ -4251,24 +4251,24 @@ dependencies = [ { name = "typing-extensions" }, ] wheels = [ - { url = "https://download-r2.pytorch.org/whl/cu126/torch-2.11.0%2Bcu126-cp311-cp311-manylinux_2_28_aarch64.whl", upload-time = "2026-03-23T14:59:31Z" }, - { url = "https://download-r2.pytorch.org/whl/cu126/torch-2.11.0%2Bcu126-cp311-cp311-manylinux_2_28_x86_64.whl", upload-time = "2026-03-23T14:59:32Z" }, - { url = "https://download-r2.pytorch.org/whl/cu126/torch-2.11.0%2Bcu126-cp311-cp311-win_amd64.whl", upload-time = "2026-03-23T14:59:35Z" }, - { url = "https://download-r2.pytorch.org/whl/cu126/torch-2.11.0%2Bcu126-cp312-cp312-manylinux_2_28_aarch64.whl", upload-time = "2026-03-23T14:59:37Z" }, - { url = "https://download-r2.pytorch.org/whl/cu126/torch-2.11.0%2Bcu126-cp312-cp312-manylinux_2_28_x86_64.whl", upload-time = "2026-03-23T14:59:37Z" }, - { url = "https://download-r2.pytorch.org/whl/cu126/torch-2.11.0%2Bcu126-cp312-cp312-win_amd64.whl", upload-time = "2026-03-23T14:59:37Z" }, - { url = "https://download-r2.pytorch.org/whl/cu126/torch-2.11.0%2Bcu126-cp313-cp313-manylinux_2_28_aarch64.whl", upload-time = "2026-03-23T14:59:44Z" }, - { url = "https://download-r2.pytorch.org/whl/cu126/torch-2.11.0%2Bcu126-cp313-cp313-manylinux_2_28_x86_64.whl", upload-time = "2026-03-23T14:59:51Z" }, - { url = "https://download-r2.pytorch.org/whl/cu126/torch-2.11.0%2Bcu126-cp313-cp313-win_amd64.whl", upload-time = "2026-03-23T14:59:51Z" }, - { url = "https://download-r2.pytorch.org/whl/cu126/torch-2.11.0%2Bcu126-cp313-cp313t-manylinux_2_28_aarch64.whl", upload-time = "2026-03-23T14:59:56Z" }, - { url = "https://download-r2.pytorch.org/whl/cu126/torch-2.11.0%2Bcu126-cp313-cp313t-manylinux_2_28_x86_64.whl", upload-time = "2026-03-23T15:00:03Z" }, - { url = "https://download-r2.pytorch.org/whl/cu126/torch-2.11.0%2Bcu126-cp313-cp313t-win_amd64.whl", upload-time = "2026-03-23T15:00:07Z" }, - { url = "https://download-r2.pytorch.org/whl/cu126/torch-2.11.0%2Bcu126-cp314-cp314-manylinux_2_28_aarch64.whl", upload-time = "2026-03-23T15:00:20Z" }, - { url = "https://download-r2.pytorch.org/whl/cu126/torch-2.11.0%2Bcu126-cp314-cp314-manylinux_2_28_x86_64.whl", upload-time = "2026-03-23T15:00:21Z" }, - { url = "https://download-r2.pytorch.org/whl/cu126/torch-2.11.0%2Bcu126-cp314-cp314-win_amd64.whl", upload-time = "2026-03-23T15:00:27Z" }, - { url = "https://download-r2.pytorch.org/whl/cu126/torch-2.11.0%2Bcu126-cp314-cp314t-manylinux_2_28_aarch64.whl", upload-time = "2026-03-23T15:00:30Z" }, - { url = "https://download-r2.pytorch.org/whl/cu126/torch-2.11.0%2Bcu126-cp314-cp314t-manylinux_2_28_x86_64.whl", upload-time = "2026-03-23T15:00:37Z" }, - { url = "https://download-r2.pytorch.org/whl/cu126/torch-2.11.0%2Bcu126-cp314-cp314t-win_amd64.whl", upload-time = "2026-03-23T15:00:39Z" }, + { url = "https://download-r2.pytorch.org/whl/cu126/torch-2.11.0%2Bcu126-cp311-cp311-manylinux_2_28_aarch64.whl", hash = "sha256:36ae7e5522b86e5d74d70f0c1a5b32eda3a0ae1635749060d0335a8856c778ea", upload-time = "2026-04-27T20:50:19Z" }, + { url = "https://download-r2.pytorch.org/whl/cu126/torch-2.11.0%2Bcu126-cp311-cp311-manylinux_2_28_x86_64.whl", hash = "sha256:a992bf8b3f2f4a8e0f2caf204ea1b1206466535eb485d4f61caf69a881e2fae7", upload-time = "2026-04-27T20:50:51Z" }, + { url = "https://download-r2.pytorch.org/whl/cu126/torch-2.11.0%2Bcu126-cp311-cp311-win_amd64.whl", hash = "sha256:ce9aeead7950ae8eb48891568f9314a9a4130996b55e797c75e0c0d2846714ae", upload-time = "2026-04-27T20:52:26Z" }, + { url = "https://download-r2.pytorch.org/whl/cu126/torch-2.11.0%2Bcu126-cp312-cp312-manylinux_2_28_aarch64.whl", hash = "sha256:dc05f49d97c04a32d139df19d77209ceca3c484083f0db1dc332f168a19ef8d8", upload-time = "2026-04-27T20:54:02Z" }, + { url = "https://download-r2.pytorch.org/whl/cu126/torch-2.11.0%2Bcu126-cp312-cp312-manylinux_2_28_x86_64.whl", hash = "sha256:9235c0c2f5032d1f8af75dd97493f3ad12aaedf81c0f3fab2d13c1d90bc54ff9", upload-time = "2026-04-27T20:54:34Z" }, + { url = "https://download-r2.pytorch.org/whl/cu126/torch-2.11.0%2Bcu126-cp312-cp312-win_amd64.whl", hash = "sha256:1b7d728ddcf84b12f64ca10ec5f6ce91dd9edcc5e62ed5f0a7260ca00f029e8b", upload-time = "2026-04-27T20:56:19Z" }, + { url = "https://download-r2.pytorch.org/whl/cu126/torch-2.11.0%2Bcu126-cp313-cp313-manylinux_2_28_aarch64.whl", hash = "sha256:3bbbfc7ed16d878aef9b91a8e78c033bf131ffb8357e23ff80911210683a31b3", upload-time = "2026-04-27T20:58:05Z" }, + { url = "https://download-r2.pytorch.org/whl/cu126/torch-2.11.0%2Bcu126-cp313-cp313-manylinux_2_28_x86_64.whl", hash = "sha256:5f6e2a5709876dd25a5fe0d4f0ab1ef33b33bb60431fcbfe768ab6468dcfe88a", upload-time = "2026-04-27T20:58:40Z" }, + { url = "https://download-r2.pytorch.org/whl/cu126/torch-2.11.0%2Bcu126-cp313-cp313-win_amd64.whl", hash = "sha256:12f7f1b73ba522576ea2050fbe2dc720a30cf681c73b8ea2d4f5480d66fad8f3", upload-time = "2026-04-27T21:00:22Z" }, + { url = "https://download-r2.pytorch.org/whl/cu126/torch-2.11.0%2Bcu126-cp313-cp313t-manylinux_2_28_aarch64.whl", hash = "sha256:7945e39fe82b85d7d018f68085ffaa6c5305a1634adb0bc792d835953b36fc9b", upload-time = "2026-04-27T21:02:11Z" }, + { url = "https://download-r2.pytorch.org/whl/cu126/torch-2.11.0%2Bcu126-cp313-cp313t-manylinux_2_28_x86_64.whl", hash = "sha256:47adcccf4de0fdb4a1550f477d919daec74f1007d15d939acfcef3080d26f9ec", upload-time = "2026-04-27T21:02:46Z" }, + { url = "https://download-r2.pytorch.org/whl/cu126/torch-2.11.0%2Bcu126-cp313-cp313t-win_amd64.whl", hash = "sha256:478ba5a57cba329f6878cf4c9bf6baa21c17c2e2b39784e0f7d3ee56b10a85a7", upload-time = "2026-04-27T21:04:33Z" }, + { url = "https://download-r2.pytorch.org/whl/cu126/torch-2.11.0%2Bcu126-cp314-cp314-manylinux_2_28_aarch64.whl", hash = "sha256:9bef6e2c04c24163bf0cafab82ee2aecdf503b1817f936e1e6f5c652231b06f1", upload-time = "2026-04-27T21:06:48Z" }, + { url = "https://download-r2.pytorch.org/whl/cu126/torch-2.11.0%2Bcu126-cp314-cp314-manylinux_2_28_x86_64.whl", hash = "sha256:09083ecca5f6f6214b021d12c971763267cd1270df9ffcbb6291c815c483b1a0", upload-time = "2026-04-27T21:07:26Z" }, + { url = "https://download-r2.pytorch.org/whl/cu126/torch-2.11.0%2Bcu126-cp314-cp314-win_amd64.whl", hash = "sha256:146a8016c5abba673ea8af7d78e45e297b23e1f53c2a6749418a7c94f5d8d9ed", upload-time = "2026-04-27T21:09:25Z" }, + { url = "https://download-r2.pytorch.org/whl/cu126/torch-2.11.0%2Bcu126-cp314-cp314t-manylinux_2_28_aarch64.whl", hash = "sha256:561ac0b9041f43c6acb927526f892cd64304eff8725536be0d86e27992588ab3", upload-time = "2026-04-27T21:11:25Z" }, + { url = "https://download-r2.pytorch.org/whl/cu126/torch-2.11.0%2Bcu126-cp314-cp314t-manylinux_2_28_x86_64.whl", hash = "sha256:e0e31c656a407c11164d2d9257e31773562ca9ad82f755851dac511d48e9fc4d", upload-time = "2026-04-27T21:12:00Z" }, + { url = "https://download-r2.pytorch.org/whl/cu126/torch-2.11.0%2Bcu126-cp314-cp314t-win_amd64.whl", hash = "sha256:3b4deea8cb01bf8f336ac8f87298f67a44835044f24deea8e0dc2d33f5d3f68d", upload-time = "2026-04-27T21:13:49Z" }, ] [[package]]