From 1707ba1d891493983919de09cdffe6b0a485903a Mon Sep 17 00:00:00 2001 From: Luca Sbardella Date: Fri, 14 Jul 2023 15:49:12 +0100 Subject: [PATCH] Fix cokpound variance --- notebooks/models/ou.md | 9 +++++++++ notebooks/models/poisson.md | 30 ++++++++++++++++++++++++++++-- notebooks/theory/levy.md | 2 +- quantflow/sp/ou.py | 9 ++++++++- quantflow/sp/poisson.py | 2 +- quantflow/utils/marginal.py | 2 +- tests/test_heston.py | 2 ++ tests/test_jump_diffusion.py | 2 +- tests/test_ou.py | 5 +++++ tests/test_poisson.py | 14 ++++++++++---- tests/utils.py | 9 +++++++++ 11 files changed, 75 insertions(+), 11 deletions(-) diff --git a/notebooks/models/ou.md b/notebooks/models/ou.md index 6279de4..45c1590 100644 --- a/notebooks/models/ou.md +++ b/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. diff --git a/notebooks/models/poisson.md b/notebooks/models/poisson.md index 8f6ab2c..44946bc 100644 --- a/notebooks/models/poisson.md +++ b/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 diff --git a/notebooks/theory/levy.md b/notebooks/theory/levy.md index ac3b4c2..d5f064e 100644 --- a/notebooks/theory/levy.md +++ b/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}} diff --git a/quantflow/sp/ou.py b/quantflow/sp/ou.py index 17f838d..3aaa20d 100755 --- a/quantflow/sp/ou.py +++ b/quantflow/sp/ou.py @@ -16,6 +16,11 @@ 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", @@ -23,7 +28,9 @@ class Vasicek(IntensityProcess): 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 diff --git a/quantflow/sp/poisson.py b/quantflow/sp/poisson.py index 09c4ccc..0341fe4 100755 --- a/quantflow/sp/poisson.py +++ b/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) diff --git a/quantflow/utils/marginal.py b/quantflow/utils/marginal.py index 571be74..6b6d463 100644 --- a/quantflow/utils/marginal.py +++ b/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() diff --git a/tests/test_heston.py b/tests/test_heston.py index de70452..0e334df 100644 --- a/tests/test_heston.py +++ b/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 diff --git a/tests/test_jump_diffusion.py b/tests/test_jump_diffusion.py index 02012fe..b39589a 100644 --- a/tests/test_jump_diffusion.py +++ b/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 diff --git a/tests/test_ou.py b/tests/test_ou.py index ee11be4..dc91de2 100644 --- a/tests/test_ou.py +++ b/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): diff --git a/tests/test_poisson.py b/tests/test_poisson.py index e401f78..032fefd 100644 --- a/tests/test_poisson.py +++ b/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,8 +36,9 @@ 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) @@ -45,7 +46,7 @@ def test_pdf(poisson: PoissonProcess) -> None: # 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 @@ -53,6 +54,11 @@ def test_sampling(poisson: PoissonProcess) -> None: 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() diff --git a/tests/utils.py b/tests/utils.py index 585b820..ee9886e 100644 --- a/tests/utils.py +++ b/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)