Fix cokpound variance

notebooks/models/ou.md

@@ -71,6 +71,15 @@ pr.sample(20, time_horizon=1, time_steps=1000).plot().update_traces(
 )
 ```
 
+```{code-cell} ipython3
+m = pr.marginal(1)
+m.mean(), m.std()
+```
+
+```{code-cell} ipython3
+m.mean_from_characteristic(), m.std_from_characteristic()
+```
+
 ## Non-gaussian OU process
 
 Non-Gaussian OU processes offer the possibility of capturing significant distributional deviations from Gaussianity and for flexible modeling of dependence structure.

notebooks/models/poisson.md

@@ -77,6 +77,12 @@ where $\Phi_{j,u}$ is the characteristic function of the jump distribution.
 
 As long as we have a closed-form solution for the characteristic function of the jump distribution, then we have a closed-form solution for the characteristic exponent of the compound Poisson process.
 
+The mean and variance of the compund Poisson is given by
+\begin{align}
+    {\mathbb E}\left[x_t\right] &= \lambda t {\mathbb E}\left[j\right]\\
+    {\mathbb Var}\left[x_t^2\right] &= \lambda t \left({\mathbb Var}\left[j\right] + {\mathbb E}\left[j\right]^2\right)
+\end{align}
+
 ### Exponential Compound Poisson Process
 
 The library includes the Exponential Poisson Process, a compound Poisson process where the jump sizes are sampled from an exponential distribution.
@@ -85,14 +91,14 @@ The library includes the Exponential Poisson Process, a compound Poisson process
 from quantflow.sp.poisson import CompoundPoissonProcess
 from quantflow.utils.distributions import Exponential
 
-pr = CompoundPoissonProcess(intensity=2, jumps=Exponential(decay=10))
+pr = CompoundPoissonProcess(intensity=1, jumps=Exponential(decay=1))
 pr
 ```
 
 ```{code-cell} ipython3
 from quantflow.utils.plot import plot_characteristic
 m = pr.marginal(1)
-plot_characteristic(m, max_frequency=100)
+plot_characteristic(m)
 ```
 
 ```{code-cell} ipython3
@@ -127,6 +133,26 @@ df = pd.DataFrame(std, index=paths.time)
 plot.plot_lines(df)
 ```
 
+## Normal Compound Poisson
+
+A compound Poisson process with a normal jump distribution
+
+```{code-cell} ipython3
+from quantflow.utils.distributions import Normal
+from quantflow.sp.poisson import CompoundPoissonProcess
+pr = CompoundPoissonProcess(intensity=10, jumps=Normal(mu=0.01, sigma=0.1))
+pr
+```
+
+```{code-cell} ipython3
+m = pr.marginal(1)
+m.mean(), m.std()
+```
+
+```{code-cell} ipython3
+m.mean_from_characteristic(), m.std_from_characteristic()
+```
+
 ## Doubly Stochastic Poisson Process
 
 The aim is to identify a stochastic process for simulating arrivals which fulfills the following properties

notebooks/theory/levy.md

@@ -28,7 +28,7 @@ Strong Markov processes. See ([Markov property](https://en.wikipedia.org/wiki/Ma
 
 ## Characteristic function
 
-The independence and stationarity of the increments of the Lévy process imply that the [characteristic function](https://en.wikipedia.org/wiki/Characteristic_function_(probability_theory)) of $x_t$ has the form
+The independence and stationarity of the increments of the Lévy process imply that the [characteristic function](./characteristic.md) of $x_t$ has the form
 
 \begin{equation}
  \Phi_{x_t, u} = {\mathbb E}\left[e^{i u x_t}\right] = e^{-\phi_{x_t, u}} = e^{-t \phi_{x_1,u}}

quantflow/sp/ou.py

@@ -16,14 +16,21 @@
 
 
 class Vasicek(IntensityProcess):
+    """Gaussian OU process, also know as the Vasiceck model.
+
+    Strictly, it is not an intensity process since it is not positive
+    """
+
     bdlp: WeinerProcess = Field(
         default_factory=WeinerProcess,
         description="Background driving Weiner process",
     )
     theta: float = Field(default=1.0, gt=0, description="Mean rate")
 
     def characteristic_exponent(self, t: FloatArrayLike, u: Vector) -> Vector:
-        raise NotImplementedError
+        mu = self.analytical_mean(t)
+        var = self.analytical_variance(t)
+        return u * (-complex(0, 1) * mu + 0.5 * var * u)
 
     def integrated_log_laplace(self, t: FloatArrayLike, u: Vector) -> Vector:
         raise NotImplementedError

quantflow/sp/poisson.py

@@ -166,4 +166,4 @@ def analytical_mean(self, t: FloatArrayLike) -> FloatArrayLike:
 
     def analytical_variance(self, t: FloatArrayLike) -> FloatArrayLike:
         """Expected variance at a time horizon"""
-        return self.intensity * t * self.jumps.variance()
+        return self.intensity * t * (self.jumps.variance() + self.jumps.mean() ** 2)

quantflow/utils/marginal.py

@@ -30,7 +30,7 @@ def support(self, points: int = 100, *, std_mult: float = 3) -> FloatArray:
     def mean(self) -> FloatArrayLike:
         """Expected value
 
-        THis should be overloaded if a more efficient way of computing the mean
+        This should be overloaded if a more efficient way of computing the mean
         """
         return self.mean_from_characteristic()
 

tests/test_heston.py

@@ -1,6 +1,7 @@
 import pytest
 
 from quantflow.sp.heston import Heston
+from tests.utils import characteristic_tests
 
 
 @pytest.fixture
@@ -12,5 +13,6 @@ def test_characteristic(heston: Heston) -> None:
     assert heston.variance_process.is_positive is True
     assert heston.characteristic(1, 0) == 1
     m = heston.marginal(1)
+    characteristic_tests(m)
     assert m.mean() == 0.0
     assert pytest.approx(m.std()) == 0.5

tests/test_jump_diffusion.py

@@ -11,7 +11,7 @@ def merton() -> Merton:
 def test_characteristic(merton: Merton) -> None:
     m = merton.marginal(1)
     assert m.mean() < 0
-    assert pytest.approx(m.std()) == 0.5
+    assert pytest.approx(m.std()) == pytest.approx(0.5, 1.0e-3)
     pdf = m.pdf_from_characteristic(128)
     assert pdf.x[0] < 0
     assert pdf.x[-1] > 0

tests/test_ou.py

@@ -2,6 +2,7 @@
 
 from quantflow.sp.bns import BNS
 from quantflow.sp.ou import GammaOU, Vasicek
+from tests.utils import characteristic_tests, analytical_tests
 
 
 @pytest.fixture
@@ -32,8 +33,12 @@ def test_sample(gamma_ou: GammaOU) -> None:
 
 def test_vasicek(vasicek: Vasicek) -> None:
     m = vasicek.marginal(10)
+    characteristic_tests(m)
+    analytical_tests(vasicek)
     assert m.mean() == 1.0
     assert m.variance() == pytest.approx(0.1)
+    assert m.mean_from_characteristic() == pytest.approx(1.0, 1e-3)
+    assert m.std_from_characteristic() == pytest.approx(m.std(), 1e-3)
 
 
 def test_bns(bns: BNS):

tests/test_poisson.py

@@ -6,7 +6,7 @@
 from quantflow.sp.dsp import DSP
 from quantflow.sp.poisson import CompoundPoissonProcess, PoissonProcess
 from quantflow.utils.distributions import Exponential
-from tests.utils import characteristic_tests
+from tests.utils import analytical_tests, characteristic_tests
 
 
 @pytest.fixture
@@ -15,7 +15,7 @@ def poisson() -> PoissonProcess:
 
 
 @pytest.fixture
-def comp() -> CompoundPoissonProcess:
+def comp() -> CompoundPoissonProcess[Exponential]:
     return CompoundPoissonProcess(intensity=2, jumps=Exponential(decay=10))
 
 
@@ -36,23 +36,29 @@ def test_characteristic(poisson: PoissonProcess) -> None:
     assert pytest.approx(m2.variance_from_characteristic(), 0.001) == 4
 
 
-def test_pdf(poisson: PoissonProcess) -> None:
+def test_poisson_pdf(poisson: PoissonProcess) -> None:
     m = poisson.marginal(1)
+    analytical_tests(poisson)
     # m.pdf(x)
     c_pdf = m.pdf_from_characteristic(32)
     np.testing.assert_almost_equal(np.linspace(0, 31, 32), c_pdf.x)
     # TODO: fix this
     # np.testing.assert_almost_equal(pdf, c_pdf.y[:10])
 
 
-def test_sampling(poisson: PoissonProcess) -> None:
+def test_poisson_sampling(poisson: PoissonProcess) -> None:
     paths = poisson.sample(1000, time_horizon=1, time_steps=1000)
     mean = paths.mean()
     assert mean[0] == 0
     std = paths.std()
     assert std[0] == 0
 
 
+def test_comp_characteristic(comp: CompoundPoissonProcess) -> None:
+    characteristic_tests(comp.marginal(1))
+    analytical_tests(comp)
+
+
 def test_dsp_sample(dsp: DSP):
     paths = dsp.sample(1000, time_horizon=1, time_steps=1000)
     mean = paths.mean()

tests/utils.py

@@ -2,6 +2,7 @@
 
 import numpy as np
 
+from quantflow.sp.base import StochasticProcess1D
 from quantflow.utils.marginal import Marginal1D
 
 
@@ -14,3 +15,11 @@ def characteristic_tests(m: Marginal1D):
     np.testing.assert_allclose(
         m.characteristic(u), cast(np.ndarray, m.characteristic(-u)).conj()
     )
+
+
+def analytical_tests(pr: StochasticProcess1D, tol: float = 1e-3):
+    t = np.linspace(0.1, 2, 20)
+    m = pr.marginal(t)
+    np.testing.assert_allclose(m.mean(), m.mean_from_characteristic(), tol)
+    np.testing.assert_allclose(m.std(), m.std_from_characteristic(), tol)
+    np.testing.assert_allclose(m.variance(), m.variance_from_characteristic(), tol)