fft test

notebooks/applications/volatility_surface.md

@@ -87,7 +87,7 @@ Moneyness is defined as
 
 ```{code-cell} ipython3
 vs.bs()
-df = vs.options_df()
+df = vs.disable_outliers(0.95).options_df()
 df
 ```
 
@@ -126,11 +126,6 @@ len(cal.options)
 cal.model
 ```
 
-```{code-cell} ipython3
-cal = cal.remove_implied_above(quantile=0.99)
-len(cal.options)
-```
-
 ```{code-cell} ipython3
 cal.fit()
 ```
@@ -140,13 +135,7 @@ pricer.model
 ```
 
 ```{code-cell} ipython3
-pricer.model.variance_process.sigma=2.8
-pricer.model.variance_process.kappa=3.96
-pricer.reset()
-```
-
-```{code-cell} ipython3
-cal.plot(index=5, max_moneyness_ttm=1)
+cal.plot(index=4, max_moneyness_ttm=1)
 ```
 
 ## Serialization

notebooks/models/heston.md

@@ -33,7 +33,7 @@ This means that the characteristic function of $y_t=x_{\tau_t}$ can be represent
 
 ```{code-cell} ipython3
 from quantflow.sp.heston import Heston
-pr = Heston.create(vol=0.6, kappa=2, sigma=1.5, rho=-0.3)
+pr = Heston.create(vol=0.6, kappa=2, sigma=1.5, rho=-0.1)
 pr
 ```
 
@@ -73,7 +73,7 @@ plot.plot_marginal_pdf(m, 128, normal=True, analytical=False, log_y=True)
 ```{code-cell} ipython3
 from quantflow.options.pricer import OptionPricer
 from quantflow.sp.heston import Heston
-pricer = OptionPricer(Heston.create(vol=0.6, kappa=2, sigma=0.8, rho=-0.1))
+pricer = OptionPricer(Heston.create(vol=0.6, kappa=2, sigma=0.8, rho=-0.2))
 pricer
 ```
 

notebooks/models/jump_diffusion.md

@@ -13,16 +13,18 @@ kernelspec:
 
 # Jump Diffusion Models
 
-The library allows to create a vast array of jump-diffusion models. The most famous one is the Merton jump diffusion model.
+The library allows to create a vast array of jump-diffusion models. The most famous one is the Merton jump-diffusion model.
+
+## Merton Model
 
 ```{code-cell} ipython3
 from quantflow.sp.jump_diffusion import Merton
 
-pr = Merton.create(diffusion_percentage=0.2, jump_intensity=50, jump_mean=0.0)
+pr = Merton.create(diffusion_percentage=0.2, jump_intensity=50, jump_skew=-0.5)
 pr
 ```
 
-## Marginal Distribution
+### Marginal Distribution
 
 ```{code-cell} ipython3
 m = pr.marginal(0.02)
@@ -35,12 +37,16 @@ from quantflow.utils import plot
 plot.plot_marginal_pdf(m, 128, normal=True, analytical=False, log_y=True)
 ```
 
-## Characteristic Function
+### Characteristic Function
 
 ```{code-cell} ipython3
 plot.plot_characteristic(m)
 ```
 
+### Option Pricing
+
+We can price options using the `OptionPricer` tooling.
+
 ```{code-cell} ipython3
 from quantflow.options.pricer import OptionPricer
 pricer = OptionPricer(pr)
@@ -54,6 +60,21 @@ for ttm in (0.05, 0.1, 0.2, 0.4, 0.6, 1):
 fig.update_layout(title="Implied black vols", height=500)
 ```
 
+This term structure of volatility demostrates one of the principal weakness of the Merton's model, and indeed of all jump diffusion models based on Lévy processes, namely the rapid flattening of the volatility surface as time-to-maturity increases.
+For very short time-to-maturities, however, the model has no problem in producing steep volatility smile and skew.
+
++++
+
+### MC paths
+
 ```{code-cell} ipython3
+pr.sample(20, time_horizon=1, time_steps=1000).plot().update_traces(line_width=0.5)
+```
+
+## Exponential Jump Diffusion
 
+This is a variation of the Mertoin model, where the jump distribution is a double exponential, one for the negative jumps and one for the positive jumps.
+
+```{code-cell} ipython3
+from 
 ```

notebooks/models/weiner.md

@@ -40,7 +40,7 @@ plot.plot_characteristic(m, n=32)
 ```{code-cell} ipython3
 from quantflow.utils import plot
 import numpy as np
-plot.plot_marginal_pdf(m, 64)
+plot.plot_marginal_pdf(m, 128)
 ```
 
 ## Test Option Pricing

notebooks/reference/glossary.md

@@ -13,6 +13,12 @@ kernelspec:
 
 # Glossary
 
+## Characteristic Function
+
+The characteristic function of a random variable $x$ is the Fourier transform of $P^x$, where $P^x$ is the distrubution measure of $x$
+\begin{equation}
+ \Phi_{x,u} = {\mathbb E}\left[e^{i u x_t}\right] = \int e^{i u x} P^x\left(dx\right)
+\end{equation}
 ## Moneyness
 
 Monenyness is used in the context of option pricing and it is defined as

quantflow/options/surface.py

@@ -395,7 +395,7 @@ def option_list(
         index: int | None = None,
         converged: bool = True,
     ) -> list[OptionPrice]:
-        "List of all option prices in the surface"
+        "List of selected option prices in the surface"
         return list(self.option_prices(select=select, index=index, converged=converged))
 
     def bs(
@@ -482,6 +482,16 @@ def as_array(
             call_put=np.array(call_put),
         )
 
+    def disable_outliers(self, quantile: float = 0.99, repeat: int = 2) -> VolSurface:
+        for _ in range(repeat):
+            option_prices = self.option_list()
+            implied_vols = [o.implied_vol for o in option_prices]
+            exclude_above = np.quantile(implied_vols, quantile)
+            for option in option_prices:
+                if option.implied_vol > exclude_above:
+                    option.converged = False
+        return self
+
     def plot(
         self,
         *,

quantflow/sp/base.py

@@ -110,7 +110,7 @@ def support(self, mean: float, std: float, points: int) -> FloatArray:
         bounds = self.domain_range()
         start = float(sigfig(lower_bound(bounds.lb, mean - std)))
         end = float(sigfig(upper_bound(bounds.ub, mean + std)))
-        return np.linspace(start, end, points)
+        return np.linspace(start, end, points + 1)
 
 
 P = TypeVar("P", bound=StochasticProcess1D)

quantflow/sp/jump_diffusion.py

@@ -36,9 +36,13 @@ def sample(self, n: int, time_horizon: float = 1, time_steps: int = 100) -> Path
         dw1 = Paths.normal_draws(n, time_horizon, time_steps)
         return self.sample_from_draws(dw1)
 
-    def sample_from_draws(self, path1: Paths, *args: Paths) -> Paths:
-        # TODO: implement
-        raise NotImplementedError
+    def sample_from_draws(self, path_w: Paths, *args: Paths) -> Paths:
+        if args:
+            path_j = args[0]
+        else:
+            path_j = self.jumps.sample(path_w.samples, path_w.t, path_w.time_steps)
+        path_w = self.diffusion.sample_from_draws(path_w)
+        return Paths(t=path_w.t, data=path_w.data + path_j.data)
 
     def analytical_mean(self, t: FloatArrayLike) -> FloatArrayLike:
         return self.diffusion.analytical_mean(t) + self.jumps.analytical_mean(t)
@@ -54,10 +58,11 @@ def create(
         vol: float = 0.5,
         diffusion_percentage: float = 0.5,
         jump_intensity: float = 100,
-        jump_mean: float = 0.0,
+        jump_skew: float = 0.0,
     ) -> Merton:
         variance = vol * vol
-        jump_std = 1
+        jump_std = 1.0
+        jump_mean = 0.0
         if diffusion_percentage > 1:
             raise ValueError("diffusion_percentage must be less than 1")
         elif diffusion_percentage < 0:
@@ -66,6 +71,7 @@ def create(
             jump_intensity = 0
         else:
             jump_std = np.sqrt(variance * (1 - diffusion_percentage) / jump_intensity)
+            jump_mean = jump_skew / jump_intensity
         return cls(
             diffusion=WeinerProcess(sigma=np.sqrt(variance * diffusion_percentage)),
             jumps=CompoundPoissonProcess[Normal](

quantflow/sp/ou.py

@@ -4,6 +4,7 @@
 
 import numpy as np
 from pydantic import Field
+from scipy.optimize import Bounds
 from scipy.stats import gamma, norm
 
 from ..utils.distributions import Exponential
@@ -45,6 +46,9 @@ def sample_from_draws(self, draws: Paths, *args: Paths) -> Paths:
             paths[t + 1, :] = x + dx
         return Paths(t=draws.t, data=paths)
 
+    def domain_range(self) -> Bounds:
+        return Bounds(-np.inf, np.inf)
+
     def analytical_mean(self, t: FloatArrayLike) -> FloatArrayLike:
         ekt = self.ekt(t)
         return self.rate * ekt + self.theta * (1 - ekt)

quantflow/utils/marginal.py

@@ -106,7 +106,7 @@ def pdf_from_characteristic(
                 simpson_rule=simpson_rule,
             )
         else:
-            x = self.support(n + 1)
+            x = self.support(n)
             min_x = float(np.min(x))
             max_x = float(np.max(x))
             delta_x = (max_x - min_x) / (len(x) - 1)

quantflow/utils/transforms.py

@@ -44,14 +44,16 @@ def default_bounds() -> Bounds:
 
 def lower_bound(b: Any, value: float) -> float:
     try:
-        return max(b[0], value)
+        v = float(b[0])
+        return value if np.isinf(v) else v
     except TypeError:
         return value
 
 
 def upper_bound(b: Any, value: float) -> float:
     try:
-        return min(b[0], value)
+        v = float(b[0])
+        return value if np.isinf(v) else v
     except TypeError:
         return value
 
@@ -94,14 +96,6 @@ def space_domain(self, delta_x: float) -> FloatArray:
             raise TransformError("Incompatible delta_x with domain bounds")
         return delta_x * grid(self.n) + b0
 
-    def delta_x_from_bounds(self) -> float | None:
-        """Return the delta_x from the domain bounds"""
-        b0 = lower_bound(self.domain_range.lb, -np.inf)
-        b1 = upper_bound(self.domain_range.ub, np.inf)
-        if np.isinf(b0) or np.isinf(b1):
-            return None
-        return (b1 - b0) / self.n
-
     def __call__(self, y: np.ndarray, delta_x: float | None = None) -> TransformResult:
         return self.fft(y) if delta_x is None else self.frft(y, delta_x)
 

tests/test_jump_diffusion.py

@@ -0,0 +1,26 @@
+import pytest
+
+from quantflow.sp.jump_diffusion import Merton
+
+
+@pytest.fixture
+def merton() -> Merton:
+    return Merton.create(diffusion_percentage=0.2, jump_skew=-0.1)
+
+
+def test_characteristic(merton: Merton) -> None:
+    m = merton.marginal(1)
+    assert m.mean() < 0
+    assert pytest.approx(m.std()) == 0.5
+    pdf = m.pdf_from_characteristic(128)
+    assert pdf.x[0] < 0
+    assert pdf.x[-1] > 0
+    assert -pdf.x[0] != pdf.x[-1]
+
+
+def test_sampling(merton: Merton) -> None:
+    paths = merton.sample(1000, time_horizon=1, time_steps=1000)
+    mean = paths.mean()
+    assert mean[0] == 0
+    std = paths.std()
+    assert std[0] == 0

tests/test_weiner.py

@@ -51,3 +51,11 @@ def test_support(weiner: WeinerProcess) -> None:
     m = weiner.marginal(0.01)
     pdf = m.pdf_from_characteristic(32)
     assert len(pdf.x) == 32
+
+
+def test_fft_v_frft(weiner: WeinerProcess) -> None:
+    m = weiner.marginal(1)
+    pdf1 = m.pdf_from_characteristic(128, max_frequency=10)
+    pdf2 = m.pdf_from_characteristic(128, use_fft=True, max_frequency=200)
+    y = np.interp(pdf1.x[10:-10], pdf2.x, pdf2.y)
+    assert np.allclose(y, pdf1.y[10:-10], 1e-2)