Fix Poisson inversion

notebooks/_toc.yml

@@ -7,6 +7,7 @@ parts:
 - caption: Theory
   chapters:
   - file: theory/levy
+  - file: theory/characteristic
   - file: theory/inversion
   - file: theory/option_pricing
 

notebooks/models/cir.md

@@ -42,7 +42,7 @@ pr.is_positive
 
 ## Marginal and moments
 
-The model has a closed-form solution for the mean, the variance, and the (marginal pdf)[https://en.wikipedia.org/wiki/Cox%E2%80%93Ingersoll%E2%80%93Ross_model].
+The model has a closed-form solution for the mean, the variance, and the [marginal pdf](https://en.wikipedia.org/wiki/Cox%E2%80%93Ingersoll%E2%80%93Ross_model).
 
 \begin{align}
 {\mathbb E}[x_t] &= x_0 e^{-\kappa t} + \theta\left(1 - e^{-\kappa t}\right) \\

notebooks/models/heston.md

@@ -127,3 +127,7 @@ std = dict(std=pr.marginal(paths.time).std(), simulated=paths.std())
 df = pd.DataFrame(std, index=paths.time)
 plot.plot_lines(df)
 ```
+
+```{code-cell} ipython3
+
+```

notebooks/models/poisson.md

@@ -82,7 +82,8 @@ As long as we have a closed-form solution for the characteristic function of the
 The library includes the Exponential Poisson Process, a compound Poisson process where the jump sizes are sampled from an exponential distribution.
 
 ```{code-cell} ipython3
-from quantflow.sp.poisson import CompoundPoissonProcess, Exponential
+from quantflow.sp.poisson import CompoundPoissonProcess
+from quantflow.utils.distributions import Exponential
 
 pr = CompoundPoissonProcess(intensity=2, jumps=Exponential(decay=10))
 pr
@@ -91,7 +92,7 @@ pr
 ```{code-cell} ipython3
 from quantflow.utils.plot import plot_characteristic
 m = pr.marginal(1)
-plot_characteristic(m)
+plot_characteristic(m, max_frequency=100)
 ```
 
 ```{code-cell} ipython3
@@ -211,6 +212,14 @@ plot_characteristic(m)
 pr.sample(10, time_horizon=10, time_steps=1000).plot().update_traces(line_width=1)
 ```
 
+```{code-cell} ipython3
+m.characteristic(2)
+```
+
+```{code-cell} ipython3
+m.characteristic(-2).conj()
+```
+
 ```{code-cell} ipython3
 
 ```

notebooks/models/weiner.md

@@ -65,5 +65,9 @@ fig.show()
 ```
 
 ```{code-cell} ipython3
-pricer.maturity(0.1, alpha=8).plot()
+pricer.maturity(0.1).plot()
+```
+
+```{code-cell} ipython3
+
 ```

notebooks/theory/characteristic.md

@@ -0,0 +1,79 @@
+---
+jupytext:
+  text_representation:
+    extension: .md
+    format_name: myst
+    format_version: 0.13
+    jupytext_version: 1.14.7
+kernelspec:
+  display_name: Python 3 (ipykernel)
+  language: python
+  name: python3
+---
+
+# Characteristic Function
+
+The library makes heavy use of characteristic function concept and therefore, it is useful to familiarize with it.
+
+## Definition
+
+The characteristic function of a random variable $x$ is the Fourier (inverse) 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}
+
+## Properties
+
+* $\Phi_{x, 0} = 1$
+* it is bounded, $\left|\Phi_{x, u}\right| \le 1$
+* it is Hermitian, $\Phi_{x, -u} = \overline{\Phi_{x, u}}$
+* it is continuous
+* characteristic function of a symmetric random variable is real-valued and even
+* moments of $x$ are given by
+\begin{equation}
+    {\mathbb E}\left[x^n\right] = i^{-n} \left.\frac{\Phi_{x, u}}{d u}\right|_{u=0}
+\end{equation}
+
+## Covolution
+
+The characteristic function is a great tool for working with linear combination of random variables.
+
+* if $x$ and $y$ are independent random variables then the characteristic function of the linear combination $a x + b y$ ($a$ and $b$ are constants) is
+
+\begin{equation}
+    \Phi_{ax+bx,u} = \Phi_{x,a u}\Phi_{y,b u}
+\end{equation}
+
+* which means, if $x$ and $y$ are independent, the characteristic function of $x+y$ is the product
+\begin{equation}
+    \Phi_{x+x,u} = \Phi_{x,u}\Phi_{y,u}
+\end{equation}
+* The characteristic function of $ax+b$ is
+\begin{equation}
+    \Phi_{ax+b,u} = e^{iub}\Phi_{x,au}
+\end{equation}
+
+## Inversion
+
+There is a one-to-one correspondence between cumulative distribution functions and characteristic functions, so it is possible to find one of these functions if we know the other.
+
+### Continuous distributions
+
+The inversion formula for these distributions is given by
+
+\begin{equation}
+    {\mathbb P}\left(x\right) = \frac{1}{2\pi}\int_{-\infty}^\infty e^{-iuk}\Phi_{x, u} du
+\end{equation}
+
+### Discrete distributions
+
+In these distributions, the random variable $x$ takes integer values. For example, the Poisson distribution is discrete.
+The inversion formula for these distributions is given by
+
+\begin{equation}
+    {\mathbb P}\left(x=k\right) = \frac{1}{2\pi}\int_{-\pi}^\pi e^{-iuk}\Phi_{x, u} du
+\end{equation}
+
+```{code-cell} ipython3
+
+```

quantflow/sp/base.py

@@ -101,9 +101,11 @@ def marginal(self, t: FloatArrayLike) -> StochasticProcess1DMarginal:
     def domain_range(self) -> Bounds:
         return default_bounds()
 
-    def max_frequency(self, std: float) -> float:
+    def frequency_range(self, std: float, max_frequency: float | None = None) -> Bounds:
         """Maximum frequency when calculating characteristic functions"""
-        return np.sqrt(40 / std / std)
+        if max_frequency is None:
+            max_frequency = np.sqrt(40 / std / std)
+        return Bounds(0, max_frequency)
 
     def support(self, mean: float, std: float, points: int) -> FloatArray:
         """Support of the process at time `t`"""
@@ -131,6 +133,10 @@ def characteristic(self, u: Vector) -> Vector:
     def domain_range(self) -> Bounds:
         return self.process.domain_range()
 
+    def frequency_range(self, max_frequency: float | None = None) -> Bounds:
+        std = float(np.min(self.std()))
+        return self.process.frequency_range(std, max_frequency=max_frequency)
+
     def pdf(self, x: FloatArrayLike) -> FloatArrayLike:
         return self.process.analytical_pdf(self.t, x)
 
@@ -149,9 +155,6 @@ def variance(self) -> FloatArrayLike:
         except NotImplementedError:
             return self.variance_from_characteristic()
 
-    def max_frequency(self) -> float:
-        return self.process.max_frequency(float(np.min(self.std())))
-
     def support(self, points: int = 100, *, std_mult: float = 4) -> FloatArray:
         return self.process.support(
             float(self.mean()), std_mult * float(self.std()), points

quantflow/sp/poisson.py

@@ -73,9 +73,9 @@ def sample_jumps(self, n: int) -> np.ndarray:
         """For a poisson process this is just a list of 1s"""
         return np.ones((n,))
 
-    def max_frequency(self, t: Vector) -> float:
-        """Maximum frequency of the process"""
-        return 2 * np.pi
+    def frequency_range(self, std: float, max_frequency: float | None = None) -> Bounds:
+        """Frequency range of the process"""
+        return Bounds(0, np.pi)
 
     def support(self, mean: float, std: float, points: int) -> FloatArray:
         """Support of the process at time `t`"""