diff --git a/Makefile b/Makefile
index cab2a5a0..e61c1864 100644
--- a/Makefile
+++ b/Makefile
@@ -31,6 +31,14 @@ book:				## Build static jupyter {book}
 	poetry run jupyter-book build notebooks --all
 
 
+.PHONY: nbconvert
+nbconvert:			## Convert notebooks to myst markdown
+	poetry run ./dev/nbconvert
+
+.PHONY: nbsync
+nbsync:				## Sync python myst notebooks to .ipynb files - needed for vs notebook development
+	poetry run ./dev/nbsync
+
 .PHONY: sphinx-config
 sphinx-config:			## Build sphinx config
 	poetry run jupyter-book config sphinx notebooks
diff --git a/dev/nbconvert b/dev/nbconvert
new file mode 100755
index 00000000..76b4453a
--- /dev/null
+++ b/dev/nbconvert
@@ -0,0 +1,7 @@
+#!/usr/bin/env bash
+set -e
+
+for file in notebooks/**/*.ipynb
+do
+  jupytext "$file" -s
+done
diff --git a/dev/nbsync b/dev/nbsync
new file mode 100755
index 00000000..d929bf95
--- /dev/null
+++ b/dev/nbsync
@@ -0,0 +1,7 @@
+#!/usr/bin/env bash
+set -e
+
+for file in notebooks/**/*.md
+do
+  jupytext "$file" -s
+done
diff --git a/notebooks/_toc.yml b/notebooks/_toc.yml
index 107ba8ef..e57c74b7 100644
--- a/notebooks/_toc.yml
+++ b/notebooks/_toc.yml
@@ -27,6 +27,7 @@ parts:
   - file: applications/overview
     sections:
     - file: applications/volatility_surface
+    - file: applications/hurst
     - file: applications/calibration
 
   - file: examples/overview
diff --git a/notebooks/applications/calibration.md b/notebooks/applications/calibration.md
index f237c7ba..9f0728f8 100644
--- a/notebooks/applications/calibration.md
+++ b/notebooks/applications/calibration.md
@@ -5,7 +5,7 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.14.7
+    jupytext_version: 1.16.6
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
@@ -26,7 +26,7 @@ Early pointers
 For calibration we use {cite:p}`ukf`.
 Lets consider the Heston model as a test case
 
-```{code-cell} ipython3
+```{code-cell}
 from quantflow.sp.heston import Heston
 
 pr = Heston.create(vol=0.6, kappa=1.3, sigma=0.8, rho=-0.6)
@@ -82,11 +82,11 @@ the state equation is given by
 X_{t+1} &= \left[\begin{matrix}\kappa\left(\theta\right) dt \\ 0\end{matrix}\right] + 
 \end{align}
 
-```{code-cell} ipython3
+```{code-cell}
 [p for p in pr.variance_process.parameters]
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 
 ```
 
@@ -106,7 +106,7 @@ x_t &= \left[\begin{matrix}\nu_t && w_t && z_t\end{matrix}\right]^T \\
 \bar{x}_t = {\mathbb E}\left[x_t\right] &= \left[\begin{matrix}\nu_t && 0 && 0\end{matrix}\right]^T
 \end{align}
 
-```{code-cell} ipython3
+```{code-cell}
 from quantflow.data.fmp import FMP
 frequency = "1min"
 async with FMP() as cli:
@@ -115,13 +115,13 @@ df = df.sort_values("date").reset_index(drop=True)
 df
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 import plotly.express as px
 fig = px.line(df, x="date", y="close", markers=True)
 fig.show()
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 import numpy as np
 from quantflow.utils.volatility import parkinson_estimator, GarchEstimator
 df["returns"] = np.log(df["close"]) - np.log(df["open"])
@@ -132,7 +132,7 @@ fig = px.line(ds["returns"], markers=True)
 fig.show()
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 import plotly.express as px
 from quantflow.utils.bins import pdf
 df = pdf(ds["returns"], num=20)
@@ -140,35 +140,35 @@ fig = px.bar(df, x="x", y="f")
 fig.show()
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 g1 = GarchEstimator.returns(ds["returns"], dt)
 g2 = GarchEstimator.pk(ds["returns"], ds["pk"], dt)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 import pandas as pd
 yf = pd.DataFrame(dict(returns=g2.y2, pk=g2.p))
 fig = px.line(yf, markers=True)
 fig.show()
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 r1 = g1.fit()
 r1
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 r2 = g2.fit()
 r2
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 sig2 = pd.DataFrame(dict(returns=np.sqrt(g2.filter(r1["params"])), pk=np.sqrt(g2.filter(r2["params"]))))
 fig = px.line(sig2, markers=False, title="Stochastic volatility")
 fig.show()
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 class HestonCalibration:
     
     def __init__(self, dt: float, initial_std = 0.5):
@@ -186,19 +186,19 @@ class HestonCalibration:
         return np.array(((1-self.kappa*self.dt, 0),(-0.5*self.dt, 0)))
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 c = HestonCalibration(dt)
 c.x0
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 c.prediction(c.x0)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 c.state_jacobian()
 ```
diff --git a/notebooks/applications/hurst.md b/notebooks/applications/hurst.md
index f234ff31..55937a4e 100644
--- a/notebooks/applications/hurst.md
+++ b/notebooks/applications/hurst.md
@@ -5,67 +5,91 @@ jupytext:
     format_name: myst
     format_version: 0.13
     jupytext_version: 1.16.6
-kernelspec:
-  display_name: Python 3 (ipykernel)
-  language: python
-  name: python3
 ---
 
-# Hurst Exponent
-
-The [Hurst exponent](https://en.wikipedia.org/wiki/Hurst_exponent) is used as a measure of long-term memory of time series. It relates to the autocorrelations of the time series, and the rate at which these decrease as the lag between pairs of values increases.
-
-It is a statistics which can be used to test if a time-series is mean reverting or it is trending.
-
-```{code-cell} ipython3
-from quantflow.sp.cir import CIR
-
-p = CIR(kappa=1, sigma=1)
-```
-
-## Study the Weiner process OHLC 
-
-```{code-cell} ipython3
+# %% [markdown]
+# # Hurst Exponent
+#
+# The [Hurst exponent](https://en.wikipedia.org/wiki/Hurst_exponent) is used as a measure of long-term memory of time series. It relates to the autocorrelations of the time series, and the rate at which these decrease as the lag between pairs of values increases.
+#
+# It is a statistics which can be used to test if a time-series is mean reverting or it is trending.
+
+# %% [markdown]
+# ## Study with the Weiner Process
+#
+# We want to construct a mechanism to estimate the Hurst exponent via OHLC data because it is widely available from data provider and easily constructed as an online signal during trading.
+#
+# In order to evaluate results against known solutions, we consider the Weiner process as generator of timeseries.
+#
+# The Weiner process is a continuous-time stochastic process named in honor of Norbert Wiener. It is often also called Brownian motion due to its historical connection with the physical model of Brownian motion of particles in water, named after the botanist Robert Brown.
+
+# %%
 from quantflow.sp.weiner import WeinerProcess
+from quantflow.utils.dates import start_of_day
 p = WeinerProcess(sigma=0.5)
-paths = p.sample(1, 1, 1000)
-df = paths.as_datetime_df().reset_index()
+paths = p.sample(1, 1, 24*60*60)
+paths.plot()
+
+# %%
+df = paths.as_datetime_df(start=start_of_day()).reset_index()
 df
-```
 
-```{code-cell} ipython3
+# %% [markdown]
+# ### Realized Variance
+#
+# At this point we estimate the standard deviation using the **realized variance** along the path (we use the **scaled** flag so that the standard deviation is scaled by the square-root of time step, in this way it removes the dependency on the time step size).
+# The value should be close to the **sigma** of the WeinerProcess defined above.
+
+# %%
+float(paths.path_std(scaled=True)[0])
+
+# %% [markdown]
+# ### Range-base Variance estimators
+#
+# We now turn our attention to range-based volatility estimators. These estimators depends on OHLC timeseries, which are widely available from data providers such as [FMP](https://site.financialmodelingprep.com/).
+# To analyze range-based variance estimators, we use he **quantflow.ta.OHLC** tool which allows to down-sample a timeserie to OHLC series and estimate variance with three different estimators
+#
+# * **Parkinson** (1980)
+# * **Garman & Klass** (1980)
+# * **Rogers & Satchell** (1991)
+#
+# See {cite:p}`molnar` for a detailed overview of the properties of range-based estimators.
+#
+# For this we build an OHLC estimator as template and use it to create OHLC estimators for different periods.
+
+# %%
+import pandas as pd
 from quantflow.ta.ohlc import OHLC
-from datetime import timedelta
 ohlc = OHLC(serie="0", period="10m", rogers_satchell_variance=True, parkinson_variance=True, garman_klass_variance=True)
-result = ohlc(df)
-result
-```
-
-```{code-cell} ipython3
-
-```
-
-# Links
-
-* [Wikipedia](https://en.wikipedia.org/wiki/Hurst_exponent)
-* [Hurst Exponent for Algorithmic Trading
-](https://robotwealth.com/demystifying-the-hurst-exponent-part-1/)
 
-```{code-cell} ipython3
+results = []
+for period in ("2m", "5m", "10m", "30m", "1h", "4h"):
+    operator = ohlc.model_copy(update=dict(period=period))
+    result = operator(df).sum()
+    results.append(dict(period=period, pk=result["0_pk"].item(), gk=result["0_gk"].item(), rs=result["0_rs"].item()))
+vdf = pd.DataFrame(results)
+vdf
+
+# %% [markdown]
+# # Links
+#
+# * [Wikipedia](https://en.wikipedia.org/wiki/Hurst_exponent)
+# * [Hurst Exponent for Algorithmic Trading
+# ](https://robotwealth.com/demystifying-the-hurst-exponent-part-1/)
+
+# %%
 import pandas as pd
 v = pd.to_timedelta(0.02, unit="d")
 v
-```
 
-```{code-cell} ipython3
+# %%
 v.to_pytimedelta()
-```
 
-```{code-cell} ipython3
+# %%
 from quantflow.utils.dates import utcnow
 pd.date_range(start=utcnow(), periods=10, freq="0.5S")
-```
 
-```{code-cell} ipython3
+# %%
+7*7+3*3
 
-```
+# %%
diff --git a/notebooks/applications/overview.md b/notebooks/applications/overview.md
index 4de72b7a..10f0567f 100644
--- a/notebooks/applications/overview.md
+++ b/notebooks/applications/overview.md
@@ -4,7 +4,7 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.14.7
+    jupytext_version: 1.16.6
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
@@ -18,6 +18,6 @@ Real-world applications of the library
 ```{tableofcontents}
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 
 ```
diff --git a/notebooks/applications/sampling.md b/notebooks/applications/sampling.md
index e2a1267c..2192c509 100644
--- a/notebooks/applications/sampling.md
+++ b/notebooks/applications/sampling.md
@@ -4,7 +4,7 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.14.7
+    jupytext_version: 1.16.6
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
@@ -15,21 +15,21 @@ kernelspec:
 
 The library use the `Paths` class for managing monte carlo paths.
 
-```{code-cell} ipython3
+```{code-cell}
 from quantflow.utils.paths import Paths
 
 nv = Paths.normal_draws(paths=1000, time_horizon=1, time_steps=1000)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 nv.var().mean()
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 nv = Paths.normal_draws(paths=1000, time_horizon=1, time_steps=1000, antithetic_variates=False)
 nv.var().mean()
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 
 ```
diff --git a/notebooks/data/fmp.md b/notebooks/data/fmp.md
index 3369a358..7a3a7d02 100644
--- a/notebooks/data/fmp.md
+++ b/notebooks/data/fmp.md
@@ -5,7 +5,7 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.16.1
+    jupytext_version: 1.16.6
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
@@ -16,45 +16,45 @@ kernelspec:
 
 The library provides a python client for the [Financial Modelling Prep API](https://site.financialmodelingprep.com/developer/docs). To use the client one needs to provide the API key aither directly to the client or via the `FMP_API_KEY` environment variable. The API offers 1 minute, 5 minutes, 15 minutes, 30 minutes, 1 hour and daily historical prices.
 
-```{code-cell} ipython3
+```{code-cell}
 from quantflow.data.fmp import FMP
 cli = FMP()
 cli.url
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 stock = "KNOS.L"
 ```
 
 ## Company Profile
 
-```{code-cell} ipython3
+```{code-cell}
 d = await cli.profile(stock)
 d
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 c = await cli.peers(stock)
 c
 ```
 
 ## Executive trading
 
-```{code-cell} ipython3
+```{code-cell}
 stock = "KNOS.L"
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 await cli.executives(stock)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 await cli.insider_trading(stock)
 ```
 
 ## News
 
-```{code-cell} ipython3
+```{code-cell}
 c = await cli.news(stock)
 c
 ```
diff --git a/notebooks/data/timeseries.md b/notebooks/data/timeseries.md
index 7729be3c..be30e752 100644
--- a/notebooks/data/timeseries.md
+++ b/notebooks/data/timeseries.md
@@ -4,7 +4,7 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.16.1
+    jupytext_version: 1.16.6
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
@@ -13,27 +13,27 @@ kernelspec:
 
 ## Timeseries
 
-```{code-cell} ipython3
+```{code-cell}
 from quantflow.data.fmp import FMP
 from quantflow.utils.plot import candlestick_plot
 cli = FMP()
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 prices = await cli.prices("ethusd", frequency="")
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 candlestick_plot(prices).update_layout(height=500)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 from quantflow.utils.df import DFutils
 
 df = DFutils(prices).with_rogers_satchel().with_parkinson()
 df
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 
 ```
diff --git a/notebooks/examples/gaussian_sampling.md b/notebooks/examples/gaussian_sampling.md
index 27d2dea3..cfb4036f 100644
--- a/notebooks/examples/gaussian_sampling.md
+++ b/notebooks/examples/gaussian_sampling.md
@@ -4,7 +4,7 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.14.7
+    jupytext_version: 1.16.6
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
@@ -15,7 +15,7 @@ kernelspec:
 
 Here we sample the gaussian OU process for different mean reversion speed and number of paths.
 
-```{code-cell} ipython3
+```{code-cell}
 from quantflow.sp.ou import Vasicek
 from quantflow.utils import plot
 import ipywidgets as widgets
@@ -51,6 +51,6 @@ fig = go.FigureWidget(data=[simulation, analytical])
 widgets.VBox([kappa, samples, fig])
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 
 ```
diff --git a/notebooks/examples/heston_vol_surface.md b/notebooks/examples/heston_vol_surface.md
index 10e9518a..2832ccf2 100644
--- a/notebooks/examples/heston_vol_surface.md
+++ b/notebooks/examples/heston_vol_surface.md
@@ -4,7 +4,7 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.14.7
+    jupytext_version: 1.16.6
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
@@ -13,7 +13,7 @@ kernelspec:
 
 # Heston Volatility Surface
 
-```{code-cell} ipython3
+```{code-cell}
 from quantflow.sp.heston import HestonJ
 from quantflow.options.pricer import OptionPricer
 
@@ -29,13 +29,13 @@ pricer = OptionPricer(model=HestonJ.create(
 pricer
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 pricer.plot3d(max_moneyness_ttm=1.5, support=31).update_layout(
     height=800,
     title="Heston volatility surface",
 )
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 
 ```
diff --git a/notebooks/examples/overview.md b/notebooks/examples/overview.md
index 3bab8fcb..738e0c17 100644
--- a/notebooks/examples/overview.md
+++ b/notebooks/examples/overview.md
@@ -4,7 +4,7 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.14.7
+    jupytext_version: 1.16.6
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
diff --git a/notebooks/examples/poisson_sampling.md b/notebooks/examples/poisson_sampling.md
index f128faa9..4dd6b211 100644
--- a/notebooks/examples/poisson_sampling.md
+++ b/notebooks/examples/poisson_sampling.md
@@ -4,7 +4,7 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.14.7
+    jupytext_version: 1.16.6
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
@@ -15,7 +15,7 @@ kernelspec:
 
 Evaluate the MC simulation for The Poisson process against the analytical PDF.
 
-```{code-cell} ipython3
+```{code-cell}
 from quantflow.sp.poisson import PoissonProcess
 from quantflow.utils import plot
 import ipywidgets as widgets
@@ -51,6 +51,6 @@ fig = go.FigureWidget(data=[simulation, analytical])
 widgets.VBox([intensity, samples, fig])
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 
 ```
diff --git a/notebooks/models/bns.md b/notebooks/models/bns.md
index 3f5774f8..cc30bf6a 100644
--- a/notebooks/models/bns.md
+++ b/notebooks/models/bns.md
@@ -4,7 +4,7 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.14.7
+    jupytext_version: 1.16.6
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
@@ -30,14 +30,14 @@ This means that the characteristic function of $y_t$ can be represented as
 
 $\phi_{w, u}$ is the characteristic exponent of $w_1$. The second equivalence is a consequence of $w$ and $\tau$ being independent, as discussed in [the time-changed Lévy](./levy.md) process section.
 
-```{code-cell} ipython3
+```{code-cell}
 from quantflow.sp.bns import BNS
 
 pr = BNS.create(vol=0.5, decay=10, kappa=10, rho=-1)
 pr
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 from quantflow.utils import plot
 m = pr.marginal(2)
 plot.plot_characteristic(m, max_frequency=10)
@@ -45,11 +45,11 @@ plot.plot_characteristic(m, max_frequency=10)
 
 ## Marginal Distribution
 
-```{code-cell} ipython3
+```{code-cell}
 m.mean(), m.std()
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 plot.plot_marginal_pdf(m, 128, normal=True, analytical=False)
 ```
 
@@ -82,11 +82,11 @@ we obtain
 
 Here we use [sympy](https://www.sympy.org/en/index.html) to derive the integral in the characteristic function.
 
-```{code-cell} ipython3
+```{code-cell}
 import sympy as sym
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 k =  sym.Symbol("k")
 iβ = sym.Symbol("iβ")
 γ = sym.Symbol("γ")
@@ -95,17 +95,17 @@ s = sym.Symbol("s")
 ϕ
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 r = sym.integrate(ϕ, s)
 sym.simplify(r)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 import numpy as np
 f = lambda x: x*np.log(x)
 f(0.001)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 
 ```
diff --git a/notebooks/models/cir.md b/notebooks/models/cir.md
index eb4511fc..fa24365b 100644
--- a/notebooks/models/cir.md
+++ b/notebooks/models/cir.md
@@ -5,7 +5,7 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.14.7
+    jupytext_version: 1.16.6
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
@@ -30,13 +30,13 @@ Importantly, the process remains positive if the Feller condition is satisfied
 
 In the code, the initial value of the process, ${\bf x}_0$, is given by the `rate` field, for example, a CIR process can be created via
 
-```{code-cell} ipython3
+```{code-cell}
 from quantflow.sp.cir import CIR
 pr = CIR(rate=1.0, kappa=2.0, sigma=1.2)
 pr
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 pr.is_positive
 ```
 
@@ -49,18 +49,18 @@ The model has a closed-form solution for the mean, the variance, and the [margin
 {\mathbb Var}[x_t] &= x_0 \frac{\sigma^2}{\kappa}\left(e^{-\kappa t} - e^{-2 \kappa t}\right) + \frac{\theta \sigma^2}{2\kappa}\left(1 - e^{-\kappa t}\right)^2 \\
 \end{align}
 
-```{code-cell} ipython3
+```{code-cell}
 m = pr.marginal(1)
 m.mean(), m.variance()
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 m.mean_from_characteristic(), m.variance_from_characteristic()
 ```
 
 The code below show the computed PDF via FRFT and the analytical formula above
 
-```{code-cell} ipython3
+```{code-cell}
 from quantflow.utils import plot
 import numpy as np
 plot.plot_marginal_pdf(m, 128, max_frequency=20)
@@ -82,7 +82,7 @@ c &= \frac{\gamma + \kappa}{2 u} \\
 d &= \frac{\gamma - \kappa}{2 u}
 \end{align}
 
-```{code-cell} ipython3
+```{code-cell}
 from quantflow.utils import plot
 m = pr.marginal(0.5)
 plot.plot_characteristic(m)
@@ -94,17 +94,17 @@ The code offers three sampling algorithms, both guarantee positiveness even if t
 
 The first sampling algorithm is the explicit Euler *full truncation* algorithm where the process is allowed to go below zero, at which point the process becomes deterministic with an upward drift of $\kappa \theta$, see {cite:p}`heston-calibration` and {cite:p}`heston-simulation` for a detailed discussion.
 
-```{code-cell} ipython3
+```{code-cell}
 from quantflow.sp.cir import CIR
 pr = CIR(rate=1.0, kappa=1.0, sigma=2.0, sample_algo="euler")
 pr
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 pr.is_positive
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 pr.sample(20, time_horizon=1, time_steps=1000).plot().update_traces(line_width=0.5)
 ```
 
@@ -112,22 +112,22 @@ The second sampling algorithm is the implicit Milstein scheme, a refinement of t
 
 The third algorithm is a fully implicit one that guarantees positiveness of the process if the Feller condition is met.
 
-```{code-cell} ipython3
+```{code-cell}
 pr = CIR(rate=1.0, kappa=1.0, sigma=0.8)
 pr
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 pr.sample(20, time_horizon=1, time_steps=1000).plot().update_traces(line_width=0.5)
 ```
 
 Sampling with a mean reversion speed 20 times larger
 
-```{code-cell} ipython3
+```{code-cell}
 pr.kappa = 20; pr
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 pr.sample(20, time_horizon=1, time_steps=1000).plot().update_traces(line_width=0.5)
 ```
 
@@ -135,7 +135,7 @@ pr.sample(20, time_horizon=1, time_steps=1000).plot().update_traces(line_width=0
 
 In this section we compare the performance of the three sampling algorithms in estimating the mean and and standard deviation.
 
-```{code-cell} ipython3
+```{code-cell}
 from quantflow.sp.cir import CIR
 
 params = dict(rate=0.8, kappa=1.5, sigma=1.2)
@@ -148,7 +148,7 @@ prs = [
 ]
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 import pandas as pd
 from quantflow.utils import plot
 from quantflow.utils.paths import Paths
@@ -164,7 +164,7 @@ df = pd.DataFrame(mean, index=draws.time)
 plot.plot_lines(df)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 std = dict(std=pr.marginal(draws.time).std())
 std.update({pr.sample_algo.name: pr.sample_from_draws(draws).std() for pr in prs})
 df = pd.DataFrame(std, index=draws.time)
diff --git a/notebooks/models/gousv.md b/notebooks/models/gousv.md
index 981b2e0b..789fc547 100644
--- a/notebooks/models/gousv.md
+++ b/notebooks/models/gousv.md
@@ -4,7 +4,7 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.14.7
+    jupytext_version: 1.16.6
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
@@ -28,7 +28,7 @@ This means that the characteristic function of $y_t=x_{\tau_t}$ can be represent
      &= e^{-a_{t,u} - b_{t,u} \nu_0}
 \end{align}
 
-```{code-cell} ipython3
+```{code-cell}
 
 ```
 
@@ -39,10 +39,10 @@ This means that the characteristic function of $y_t=x_{\tau_t}$ can be represent
     b_t &= \frac{1 - e^{-\kappa t}}{\kappa} \\
 \end{align}
 
-```{code-cell} ipython3
+```{code-cell}
 
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 
 ```
diff --git a/notebooks/models/heston.md b/notebooks/models/heston.md
index b2c926fb..f012aa96 100644
--- a/notebooks/models/heston.md
+++ b/notebooks/models/heston.md
@@ -5,7 +5,7 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.14.7
+    jupytext_version: 1.16.6
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
@@ -31,20 +31,20 @@ This means that the characteristic function of $y_t=x_{\tau_t}$ can be represent
      &= e^{-a_{t,u} - b_{t,u} \nu_0}
 \end{align}
 
-```{code-cell} ipython3
+```{code-cell}
 from quantflow.sp.heston import Heston
 pr = Heston.create(vol=0.6, kappa=2, sigma=1.5, rho=-0.1)
 pr
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 # check that the variance CIR process is positive
 pr.variance_process.is_positive, pr.variance_process.marginal(1).std()
 ```
 
 ## Characteristic Function
 
-```{code-cell} ipython3
+```{code-cell}
 from quantflow.utils import plot
 m = pr.marginal(0.1)
 plot.plot_characteristic(m)
@@ -58,26 +58,26 @@ The immaginary part of the characteristic function is given by the correlation c
 
 Here we compare the marginal distribution at a time in the future $t=1$ with a normal distribution with the same standard deviation.
 
-```{code-cell} ipython3
+```{code-cell}
 plot.plot_marginal_pdf(m, 128, normal=True, analytical=False)
 ```
 
 Using log scale on the y axis highlighs the probability on the tails much better
 
-```{code-cell} ipython3
+```{code-cell}
 plot.plot_marginal_pdf(m, 128, normal=True, analytical=False, log_y=True)
 ```
 
 ## Option pricing
 
-```{code-cell} ipython3
+```{code-cell}
 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.2))
 pricer
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 import plotly.express as px
 import plotly.graph_objects as go
 from quantflow.options.bs import black_call
@@ -89,7 +89,7 @@ fig.add_trace(go.Scatter(x=r.moneyness_ttm, y=b.time_value, name=b.name, line=di
 fig.show()
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 fig = None
 for ttm in (0.05, 0.1, 0.2, 0.4, 0.6, 1):
     fig = pricer.maturity(ttm).plot(fig=fig, name=f"t={ttm}")
@@ -102,17 +102,17 @@ The simulation of the Heston model is heavily dependent on the simulation of the
 
 The code implements algorithms from {cite:p}heston-simulation
 
-```{code-cell} ipython3
+```{code-cell}
 from quantflow.sp.heston import Heston
 pr = Heston.create(vol=0.6, kappa=2, sigma=0.8, rho=-0.4)
 pr
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 pr.sample(20, time_horizon=1, time_steps=1000).plot().update_traces(line_width=0.5)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 import pandas as pd
 from quantflow.utils import plot
 
@@ -122,12 +122,12 @@ df = pd.DataFrame(mean, index=paths.time)
 plot.plot_lines(df)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 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
+```{code-cell}
 
 ```
diff --git a/notebooks/models/heston_jumps.md b/notebooks/models/heston_jumps.md
index 0c1a021a..10012a9c 100644
--- a/notebooks/models/heston_jumps.md
+++ b/notebooks/models/heston_jumps.md
@@ -4,7 +4,7 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.14.7
+    jupytext_version: 1.16.6
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
@@ -32,7 +32,7 @@ where $j_t$ is a double exponential Compound Poisson process which adds three ad
 
 The jump process is independent of the other Brownian motions.
 
-```{code-cell} ipython3
+```{code-cell}
 from quantflow.sp.heston import HestonJ
 pr = HestonJ.create(
     vol=0.6,
@@ -46,25 +46,25 @@ pr = HestonJ.create(
 pr
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 from quantflow.utils import plot
 plot.plot_marginal_pdf(pr.marginal(0.1), 128, normal=True, analytical=False)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 from quantflow.options.pricer import OptionPricer
 from quantflow.sp.heston import Heston
 pricer = OptionPricer(pr)
 pricer
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 fig = None
 for ttm in (0.05, 0.1, 0.2, 0.4, 0.6, 1):
     fig = pricer.maturity(ttm).plot(fig=fig, name=f"t={ttm}")
 fig.update_layout(title="Implied black vols", height=500)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 
 ```
diff --git a/notebooks/models/jump_diffusion.md b/notebooks/models/jump_diffusion.md
index f06079da..23dd3fd8 100644
--- a/notebooks/models/jump_diffusion.md
+++ b/notebooks/models/jump_diffusion.md
@@ -4,7 +4,7 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.14.7
+    jupytext_version: 1.16.6
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
@@ -17,7 +17,7 @@ The library allows to create a vast array of jump-diffusion models. The most fam
 
 ## Merton Model
 
-```{code-cell} ipython3
+```{code-cell}
 from quantflow.sp.jump_diffusion import Merton
 
 pr = Merton.create(diffusion_percentage=0.2, jump_intensity=50, jump_skew=-0.5)
@@ -26,12 +26,12 @@ pr
 
 ### Marginal Distribution
 
-```{code-cell} ipython3
+```{code-cell}
 m = pr.marginal(0.02)
 m.std(), m.std_from_characteristic()
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 from quantflow.utils import plot
 
 plot.plot_marginal_pdf(m, 128, normal=True, analytical=False, log_y=True)
@@ -39,7 +39,7 @@ plot.plot_marginal_pdf(m, 128, normal=True, analytical=False, log_y=True)
 
 ### Characteristic Function
 
-```{code-cell} ipython3
+```{code-cell}
 plot.plot_characteristic(m)
 ```
 
@@ -47,13 +47,13 @@ plot.plot_characteristic(m)
 
 We can price options using the `OptionPricer` tooling.
 
-```{code-cell} ipython3
+```{code-cell}
 from quantflow.options.pricer import OptionPricer
 pricer = OptionPricer(pr)
 pricer
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 fig = None
 for ttm in (0.05, 0.1, 0.2, 0.4, 0.6, 1):
     fig = pricer.maturity(ttm).plot(fig=fig, name=f"t={ttm}")
@@ -67,7 +67,7 @@ For very short time-to-maturities, however, the model has no problem in producin
 
 ### MC paths
 
-```{code-cell} ipython3
+```{code-cell}
 pr.sample(20, time_horizon=1, time_steps=1000).plot().update_traces(line_width=0.5)
 ```
 
@@ -75,6 +75,6 @@ pr.sample(20, time_horizon=1, time_steps=1000).plot().update_traces(line_width=0
 
 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
+```{code-cell}
 from 
 ```
diff --git a/notebooks/models/ou.md b/notebooks/models/ou.md
index 45c15901..3a2ce252 100644
--- a/notebooks/models/ou.md
+++ b/notebooks/models/ou.md
@@ -4,7 +4,7 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.14.7
+    jupytext_version: 1.16.6
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
@@ -57,13 +57,13 @@ which means the process admits a stationary probability distribution equal to
     x_t \sim N\left(\theta, \frac{\sigma^2}{2\kappa}\right)\ \ t\rightarrow\infty
 \end{equation}
 
-```{code-cell} ipython3
+```{code-cell}
 from quantflow.sp.ou import Vasicek
 pr = Vasicek()
 pr
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 pr.sample(20, time_horizon=1, time_steps=1000).plot().update_traces(
     line_width=0.5
 ).update_layout(
@@ -71,12 +71,12 @@ pr.sample(20, time_horizon=1, time_steps=1000).plot().update_traces(
 )
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 m = pr.marginal(1)
 m.mean(), m.std()
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 m.mean_from_characteristic(), m.std_from_characteristic()
 ```
 
@@ -125,7 +125,7 @@ The library provides an implementation of the non-gaussian OU process in the for
 
 In this case, the BDLP is an exponential compound Poisson process with Lévy density $\lambda\beta e^{-\beta x}$, in other words, the [exponential compound Poisson](./poisson.md) process with intensity $\lambda$ and decay $\beta$.
 
-```{code-cell} ipython3
+```{code-cell}
 from quantflow.sp.ou import GammaOU
 
 pr = GammaOU.create(decay=10, kappa=5)
@@ -140,11 +140,11 @@ The charatecristic exponent of the $\Gamma$-OU process is given by, see {cite:p}
     \phi_{u, t} = -x_{0} i u e^{-\kappa t} - \lambda\ln\left(\frac{\beta-iue^{-\kappa t}}{\beta -iu}\right)
 \end{equation}
 
-```{code-cell} ipython3
+```{code-cell}
 pr.marginal(1).mean(), pr.marginal(1).std()
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 import numpy as np
 from quantflow.utils import plot
 
@@ -152,14 +152,14 @@ m = pr.marginal(5)
 plot.plot_marginal_pdf(m, 128)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 from quantflow.utils.plot import plot_characteristic
 plot_characteristic(m)
 ```
 
 ### Sampling Gamma OU
 
-```{code-cell} ipython3
+```{code-cell}
 from quantflow.sp.ou import GammaOU
 pr = GammaOU.create(decay=10, kappa=5)
 
@@ -170,7 +170,7 @@ pr.sample(50, time_horizon=1, time_steps=1000).plot().update_traces(line_width=0
 
 Test the simulated meand and stadard deviation against the values from the invariant gamma distribution.
 
-```{code-cell} ipython3
+```{code-cell}
 import pandas as pd
 from quantflow.utils import plot
 
@@ -180,7 +180,7 @@ df = pd.DataFrame(mean, index=paths.time)
 plot.plot_lines(df)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 std = dict(std=pr.marginal(paths.time).std(), simulated=paths.std())
 df = pd.DataFrame(std, index=paths.time)
 plot.plot_lines(df)
@@ -198,6 +198,6 @@ The integration of the OU process can be achieved by multiplying both sides of t
     x_t &= x_0 e^{-\kappa t} + \int_0^t e^{-\kappa\left(t - s\right)} d z_s
 \end{align}
 
-```{code-cell} ipython3
+```{code-cell}
 
 ```
diff --git a/notebooks/models/overview.md b/notebooks/models/overview.md
index 541ee036..704d05f3 100644
--- a/notebooks/models/overview.md
+++ b/notebooks/models/overview.md
@@ -4,7 +4,7 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.14.7
+    jupytext_version: 1.16.6
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
diff --git a/notebooks/models/poisson.md b/notebooks/models/poisson.md
index c5c71d91..3e0acd3a 100644
--- a/notebooks/models/poisson.md
+++ b/notebooks/models/poisson.md
@@ -5,7 +5,7 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.14.7
+    jupytext_version: 1.16.6
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
@@ -31,25 +31,25 @@ The characteristic exponent is given by
 \phi_{N_t, u} = t \lambda \left(1 - e^{iu}\right)
 \end{equation}
 
-```{code-cell} ipython3
+```{code-cell}
 from quantflow.sp.poisson import PoissonProcess
 pr = PoissonProcess(intensity=1)
 pr
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 m = pr.marginal(0.1)
 import numpy as np
 cdf = m.cdf(np.arange(5))
 cdf
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 n = 128*8
 m.cdf_from_characteristic(5, frequency_n=n).y
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 cdf1 = m.cdf_from_characteristic(5, frequency_n=n).y
 cdf2 = m.cdf_from_characteristic(5, frequency_n=n, simpson_rule=False).y
 10000*np.max(np.abs(cdf-cdf1)), 10000*np.max(np.abs(cdf-cdf2))
@@ -57,7 +57,7 @@ cdf2 = m.cdf_from_characteristic(5, frequency_n=n, simpson_rule=False).y
 
 ### Marginal
 
-```{code-cell} ipython3
+```{code-cell}
 import numpy as np
 from quantflow.utils import plot
 
@@ -65,14 +65,14 @@ m = pr.marginal(1)
 plot.plot_marginal_pdf(m, frequency_n=128*8)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 from quantflow.utils.plot import plot_characteristic
 plot_characteristic(m)
 ```
 
 ### Sampling Poisson
 
-```{code-cell} ipython3
+```{code-cell}
 p = pr.sample(10, time_horizon=10, time_steps=1000)
 p.plot().update_traces(line_width=1)
 ```
@@ -105,7 +105,7 @@ The mean and variance of the compund Poisson is given by
 
 The library includes the Exponential Poisson Process, a compound Poisson process where the jump sizes are sampled from an exponential distribution.
 
-```{code-cell} ipython3
+```{code-cell}
 from quantflow.sp.poisson import CompoundPoissonProcess
 from quantflow.utils.distributions import Exponential
 
@@ -113,21 +113,21 @@ pr = CompoundPoissonProcess(intensity=1, jumps=Exponential(decay=1))
 pr
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 from quantflow.utils.plot import plot_characteristic
 m = pr.marginal(1)
 plot_characteristic(m)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 m.mean(), m.mean_from_characteristic()
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 m.variance(), m.variance_from_characteristic()
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 pr.sample(10, time_horizon=10, time_steps=1000).plot().update_traces(line_width=1)
 ```
 
@@ -135,7 +135,7 @@ pr.sample(10, time_horizon=10, time_steps=1000).plot().update_traces(line_width=
 
 Here we test the simulated mean and standard deviation against the analytical values.
 
-```{code-cell} ipython3
+```{code-cell}
 import pandas as pd
 from quantflow.utils import plot
 
@@ -145,7 +145,7 @@ df = pd.DataFrame(mean, index=paths.time)
 plot.plot_lines(df)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 std = dict(std=pr.marginal(paths.time).std(), simulated=paths.std())
 df = pd.DataFrame(std, index=paths.time)
 plot.plot_lines(df)
@@ -155,19 +155,19 @@ plot.plot_lines(df)
 
 A compound Poisson process with a normal jump distribution
 
-```{code-cell} ipython3
+```{code-cell}
 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
+```{code-cell}
 m = pr.marginal(1)
 m.mean(), m.std()
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 m.mean_from_characteristic(), m.std_from_characteristic()
 ```
 
@@ -222,14 +222,14 @@ The intensity function of a DSPP is given by:
 {\mathbb P}\left(N_T - N_t = n\right) = {\mathbb E}_t\left[e^{-\Lambda_{t,T}} \frac{\Lambda_{t, T}^n}{n!}\right] = \frac{1}{n!}
 \end{equation}
 
-```{code-cell} ipython3
+```{code-cell}
 from quantflow.sp.dsp import DSP, PoissonProcess, CIR
 pr = DSP(intensity=CIR(sigma=2, kappa=1), poisson=PoissonProcess(intensity=2))
 pr2 = DSP(intensity=CIR(rate=2, sigma=4, kappa=2, theta=2), poisson=PoissonProcess(intensity=1))
 pr, pr2
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 import numpy as np
 from quantflow.utils import plot
 import plotly.graph_objects as go
@@ -242,28 +242,28 @@ plot.plot_marginal_pdf(pr2.marginal(1), n, analytical=False, fig=fig, marker_col
 fig.add_trace(go.Scatter(x=pdf.x, y=pr.poisson.marginal(1).pdf(pdf.x), name="Poisson", mode="markers", marker_color="blue"))
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 pr.marginal(1).mean(), pr.marginal(1).variance()
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 pr2.marginal(1).mean(), pr2.marginal(1).variance()
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 from quantflow.utils.plot import plot_characteristic
 m = pr.marginal(2)
 plot_characteristic(m)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 pr.sample(10, time_horizon=10, time_steps=1000).plot().update_traces(line_width=1)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 m.characteristic(2)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 m.characteristic(-2).conj()
 ```
diff --git a/notebooks/models/weiner.md b/notebooks/models/weiner.md
index 75de1652..0d496448 100644
--- a/notebooks/models/weiner.md
+++ b/notebooks/models/weiner.md
@@ -4,7 +4,7 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.14.7
+    jupytext_version: 1.16.6
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
@@ -23,21 +23,21 @@ The characteristic exponent of $w$ is
     \phi_{w, u} = \frac{\sigma^2 u^2}{2}
 \end{equation}
 
-```{code-cell} ipython3
+```{code-cell}
 from quantflow.sp.weiner import WeinerProcess
 
 pr = WeinerProcess(sigma=0.5)
 pr
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 from quantflow.utils import plot
 # create the marginal at time in the future
 m = pr.marginal(1)
 plot.plot_characteristic(m, n=32)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 from quantflow.utils import plot
 import numpy as np
 plot.plot_marginal_pdf(m, 128)
@@ -45,14 +45,14 @@ plot.plot_marginal_pdf(m, 128)
 
 ## Test Option Pricing
 
-```{code-cell} ipython3
+```{code-cell}
 from quantflow.options.pricer import OptionPricer
 from quantflow.sp.weiner import WeinerProcess
 pricer = OptionPricer(WeinerProcess(sigma=0.2))
 pricer
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 import plotly.express as px
 import plotly.graph_objects as go
 from quantflow.options.bs import black_call
@@ -64,10 +64,10 @@ fig.add_trace(go.Scatter(x=b.moneyness_ttm, y=b.time_value, name=b.name, line=di
 fig.show()
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 pricer.maturity(0.1).plot()
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 
 ```
diff --git a/notebooks/reference/biblio.md b/notebooks/reference/biblio.md
index 75560014..7e54e04f 100644
--- a/notebooks/reference/biblio.md
+++ b/notebooks/reference/biblio.md
@@ -5,7 +5,7 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.13.8
+    jupytext_version: 1.16.6
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
diff --git a/notebooks/reference/contributing.md b/notebooks/reference/contributing.md
index 7d7ab288..f1c8ba8c 100644
--- a/notebooks/reference/contributing.md
+++ b/notebooks/reference/contributing.md
@@ -5,7 +5,7 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.14.7
+    jupytext_version: 1.16.6
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
diff --git a/notebooks/reference/glossary.md b/notebooks/reference/glossary.md
index 0f1bb21d..7bcc9013 100644
--- a/notebooks/reference/glossary.md
+++ b/notebooks/reference/glossary.md
@@ -4,7 +4,7 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.14.7
+    jupytext_version: 1.16.6
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
@@ -29,6 +29,6 @@ Monenyness is used in the context of option pricing and it is defined as
 
 where $K$ is the strike and $F$ is the Forward price. A positive value implies strikes above the forward, which means put options are in the money and call options are out of the money.
 
-```{code-cell} ipython3
+```{code-cell}
 
 ```
diff --git a/notebooks/reference/references.bib b/notebooks/reference/references.bib
index 9ebc0236..0aec0d75 100644
--- a/notebooks/reference/references.bib
+++ b/notebooks/reference/references.bib
@@ -130,3 +130,10 @@ @article{dspp
     year={2017},
 }
 
+@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},
+}
diff --git a/notebooks/theory/characteristic.md b/notebooks/theory/characteristic.md
index 5c73d781..37207089 100644
--- a/notebooks/theory/characteristic.md
+++ b/notebooks/theory/characteristic.md
@@ -4,7 +4,7 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.14.7
+    jupytext_version: 1.16.6
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
@@ -74,6 +74,6 @@ The inversion formula for these distributions is given by
     {\mathbb P}\left(x=k\right) = \frac{1}{2\pi}\int_{-\pi}^\pi e^{-iuk}\Phi_{x, u} du
 \end{equation}
 
-```{code-cell} ipython3
+```{code-cell}
 
 ```
diff --git a/notebooks/theory/inversion.md b/notebooks/theory/inversion.md
index 1e4cfe14..e147cee1 100644
--- a/notebooks/theory/inversion.md
+++ b/notebooks/theory/inversion.md
@@ -5,7 +5,7 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.14.7
+    jupytext_version: 1.16.6
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
@@ -72,26 +72,26 @@ which means $\delta_u$ and $\delta_x$ cannot be chosen indipendently.
 
 As an example, let us invert the characteristic function of the Weiber process, which yields the standard distribution.
 
-```{code-cell} ipython3
+```{code-cell}
 from quantflow.sp.weiner import WeinerProcess
 p = WeinerProcess(sigma=0.5)
 m = p.marginal(0.2)
 m.std()
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 from quantflow.utils import plot
 
 plot.plot_characteristic(m)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 from quantflow.utils import plot
 import numpy as np
 plot.plot_marginal_pdf(m, 128, use_fft=True, max_frequency=20)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 plot.plot_marginal_pdf(m, 128*8, use_fft=True, max_frequency=8*20)
 ```
 
@@ -113,7 +113,7 @@ z &= \left(\left[e^{i j^2 \zeta/2}\right]_{j=0}^{N-1}, \left[e^{i\left(N-j\right
 
 We can now reduce the number of points needed for the discretization and achieve higher accuracy by properly selecting the domain discretization independently.
 
-```{code-cell} ipython3
+```{code-cell}
 plot.plot_marginal_pdf(m, 128)
 ```
 
diff --git a/notebooks/theory/levy.md b/notebooks/theory/levy.md
index d5f064e7..a45664af 100644
--- a/notebooks/theory/levy.md
+++ b/notebooks/theory/levy.md
@@ -6,7 +6,7 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.14.7
+    jupytext_version: 1.16.6
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
@@ -116,6 +116,6 @@ When the intensity process is affine, the Laplace transform takes the following
 
 where coefficients $a$ and $b$ satisfy Riccati ODEs, which can be solved numerically and, in some cases, analytically.
 
-```{code-cell} ipython3
+```{code-cell}
 
 ```
diff --git a/notebooks/theory/option_pricing.md b/notebooks/theory/option_pricing.md
index 5898c6a1..b4b32997 100644
--- a/notebooks/theory/option_pricing.md
+++ b/notebooks/theory/option_pricing.md
@@ -5,7 +5,7 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.14.7
+    jupytext_version: 1.16.6
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
@@ -48,7 +48,7 @@ This is the very famous convexity correction which appears in all diffusion driv
     c_t = \frac{\sigma^2 t}{2}
 \end{equation}
 
-```{code-cell} ipython3
+```{code-cell}
 from quantflow.sp.weiner import WeinerProcess
 pr = WeinerProcess(sigma=0.5)
 -pr.characteristic_exponent(1, complex(0,-1))
@@ -56,7 +56,7 @@ pr = WeinerProcess(sigma=0.5)
 
 which is the same as
 
-```{code-cell} ipython3
+```{code-cell}
 pr.convexity_correction(1)
 ```
 
@@ -109,7 +109,7 @@ is provided by $\Psi_{-i\alpha}$ being finite, which means the characteristic fu
 
 Here we illustrate how to use the characteristic function integration with the classical [Weiner process](https://en.wikipedia.org/wiki/Wiener_process).
 
-```{code-cell} ipython3
+```{code-cell}
 from quantflow.sp.weiner import WeinerProcess
 ttm=1
 p = WeinerProcess(sigma=0.5)
@@ -119,7 +119,7 @@ m = p.marginal(ttm)
 m.std()
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 import plotly.express as px
 import plotly.graph_objects as go
 from quantflow.options.bs import black_call
@@ -132,6 +132,6 @@ fig.add_trace(go.Scatter(x=r.x, y=b, name="analytical", line=dict()))
 fig.show()
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 
 ```
diff --git a/notebooks/theory/overview.md b/notebooks/theory/overview.md
index 1d0e1b06..b75da9f9 100644
--- a/notebooks/theory/overview.md
+++ b/notebooks/theory/overview.md
@@ -4,7 +4,7 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.14.7
+    jupytext_version: 1.16.6
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
diff --git a/poetry.lock b/poetry.lock
index 3a46e9da..a6327285 100644
--- a/poetry.lock
+++ b/poetry.lock
@@ -4798,4 +4798,4 @@ data = ["aio-fluid"]
 [metadata]
 lock-version = "2.0"
 python-versions = ">=3.11"
-content-hash = "945f85a4526b83f35e916bb16c00d24c25e34667b1e606ddf906ef348c73b2e8"
+content-hash = "1e04c99ee3e0472ff9fde2a84cb1d44554699bf5837b7a1f83ff3eb1d53db145"
diff --git a/pyproject.toml b/pyproject.toml
index dcd58d5e..2a4305ee 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -43,7 +43,6 @@ optional = true
 
 [tool.poetry.group.book.dependencies]
 jupyter-book = "^1.0.0"
-nbconvert = "^7.16.3"
 jupytext = "^1.13.8"
 plotly = "^5.20.0"
 jupyterlab = "^4.0.2"
diff --git a/quantflow/data/vault.py b/quantflow/data/vault.py
index d54f90b1..0178aa4f 100644
--- a/quantflow/data/vault.py
+++ b/quantflow/data/vault.py
@@ -2,6 +2,7 @@
 
 
 class Vault:
+    """Keeps key-value pairs in a file."""
 
     def __init__(self, path: str | Path) -> None:
         self.path = Path(path)
@@ -17,22 +18,27 @@ def load(self) -> dict[str, str]:
         return data
 
     def add(self, key: str, value: str) -> None:
+        """Add a key-value pair to the vault."""
         self.data[key] = value
         self.save()
 
     def delete(self, key: str) -> bool:
+        """Delete a key-value pair from the vault."""
         if self.data.pop(key, None) is not None:
             self.save()
             return True
         return False
 
     def get(self, key: str) -> str | None:
+        """Get the value of a key if available otherwise None."""
         return self.data.get(key)
 
     def keys(self) -> list[str]:
+        """Get the keys in the vault."""
         return sorted(self.data)
 
     def save(self) -> None:
+        """Save the data to the file."""
         with open(self.path, "w") as file:
             for key in sorted(self.data):
                 value = self.data[key]
diff --git a/quantflow/ta/base.py b/quantflow/ta/base.py
index 47ea8b8a..a49bc568 100644
--- a/quantflow/ta/base.py
+++ b/quantflow/ta/base.py
@@ -6,7 +6,9 @@
 DataFrame: TypeAlias = pl.DataFrame | pd.DataFrame
 
 
-def to_polars(df: DataFrame) -> pl.DataFrame:
+def to_polars(df: DataFrame, *, copy: bool = False) -> pl.DataFrame:
     if isinstance(df, pd.DataFrame):
         return pl.DataFrame(df)
+    elif copy:
+        return df.clone()
     return df
diff --git a/quantflow/ta/ohlc.py b/quantflow/ta/ohlc.py
index 86291d91..ab8c1274 100644
--- a/quantflow/ta/ohlc.py
+++ b/quantflow/ta/ohlc.py
@@ -1,14 +1,13 @@
-from dataclasses import dataclass
 from datetime import timedelta
 
 import numpy as np
 import polars as pl
+from pydantic import BaseModel
 
 from .base import DataFrame, to_polars
 
 
-@dataclass
-class OHLC:
+class OHLC(BaseModel):
     """Aggregates OHLC data over a given period and serie
 
     Optionally calculates the range-based variance estimators for the serie.
@@ -50,7 +49,7 @@ def close_col(self) -> pl.Expr:
     def __call__(self, df: DataFrame) -> pl.DataFrame:
         """Returns a dataframe with OHLC data sampled over the given period"""
         result = (
-            to_polars(df)
+            to_polars(df, copy=True)
             .group_by_dynamic(self.index_column, every=self.period)
             .agg(
                 pl.col(self.serie).first().alias(f"{self.serie}_open"),
diff --git a/quantflow/utils/dates.py b/quantflow/utils/dates.py
index de012576..63829c93 100644
--- a/quantflow/utils/dates.py
+++ b/quantflow/utils/dates.py
@@ -5,7 +5,20 @@ def utcnow() -> datetime:
     return datetime.now(timezone.utc)
 
 
+def as_utc(dt: date | None = None) -> datetime:
+    if dt is None:
+        return utcnow()
+    elif isinstance(dt, datetime):
+        return dt.astimezone(timezone.utc)
+    else:
+        return datetime(dt.year, dt.month, dt.day, tzinfo=timezone.utc)
+
+
 def isoformat(date: str | date) -> str:
     if isinstance(date, str):
         return date
     return date.isoformat()
+
+
+def start_of_day(dt: date | None = None) -> datetime:
+    return as_utc(dt).replace(hour=0, minute=0, second=0, microsecond=0)
diff --git a/quantflow/utils/paths.py b/quantflow/utils/paths.py
index 42b6dc52..ebf18a62 100755
--- a/quantflow/utils/paths.py
+++ b/quantflow/utils/paths.py
@@ -19,10 +19,13 @@ class Paths(BaseModel, arbitrary_types_allowed=True):
     """Paths of a stochastic process"""
 
     t: float = Field(description="time horizon")
+    """Time horizon - the unit of time is not specified"""
     data: FloatArray = Field(description="paths")
+    """Paths of the stochastic process"""
 
     @property
     def dt(self) -> float:
+        """Time step"""
         return self.t / self.time_steps
 
     @property
@@ -64,17 +67,42 @@ def dates(
         return pd.date_range(start=start, end=end, periods=self.time_steps + 1)
 
     def mean(self) -> FloatArray:
-        """Mean of paths"""
+        """Paths cross-section mean"""
         return np.mean(self.data, axis=1)
 
     def std(self) -> FloatArray:
-        """Standard deviation of paths"""
+        """Paths cross-section standard deviation"""
         return np.std(self.data, axis=1)
 
     def var(self) -> FloatArray:
-        """Variance of paths"""
+        """Paths cross-section variance"""
         return np.var(self.data, axis=1)
 
+    def paths_mean(self, *, scaled: bool = False) -> FloatArray:
+        """mean for each path
+
+        If scaled is True, the mean is scaled by the time step
+        """
+        scale = self.dt if scaled else 1.0
+        return np.mean(self.data, axis=0) / scale
+
+    def paths_std(self, *, scaled: bool = False) -> FloatArray:
+        """standard deviation for each path
+
+        If scaled is True, the standard deviation is scaled by the square
+        root of the time step
+        """
+        scale = np.sqrt(self.dt) if scaled else 1.0
+        return np.std(np.diff(self.data, axis=0), axis=0) / scale
+
+    def paths_var(self, *, scaled: bool = False) -> FloatArray:
+        """variance for each path
+
+        If scaled is True, the variance is scaled by the time step
+        """
+        scale = self.dt if scaled else 1.0
+        return np.var(np.diff(self.data, axis=0), axis=0) / scale
+
     def as_datetime_df(
         self, *, start: datetime | None = None, unit: str = "d"
     ) -> pd.DataFrame:
diff --git a/quantflow/utils/volatility.py b/quantflow/utils/volatility.py
deleted file mode 100644
index deb945e5..00000000
--- a/quantflow/utils/volatility.py
+++ /dev/null
@@ -1,71 +0,0 @@
-from typing import Any
-
-import numpy as np
-from scipy.optimize import minimize
-
-from .types import Vector
-
-
-def parkinson_estimator(high: Vector, low: Vector) -> Vector:
-    """Parkinson volatility estimator"""
-    return np.power(np.log(high) - np.log(low), 2) / np.log(2) / 4
-
-
-def garman_klass_estimator(high: Vector, low: Vector, close: Vector) -> Vector:
-    """Garman-Klass volatility estimator"""
-    return 0.5 * np.power(np.log(high) - np.log(low), 2) - (
-        2 * np.log(2) - 1
-    ) * np.power(np.log(close), 2)
-
-
-def akaike_information(p: np.ndarray, log_like: float) -> float:
-    return 2 * p.shape[0] - 2 * log_like
-
-
-class GarchEstimator:
-    """GARCH 1,1 volatility estimator"""
-
-    def __init__(self, y2: np.ndarray, p: np.ndarray):
-        self.y2 = y2
-        self.p = p
-
-    @property
-    def n(self) -> int:
-        return self.y2.shape[0]
-
-    @classmethod
-    def returns(cls, y: np.ndarray, dt: float = 1.0) -> "GarchEstimator":
-        y = np.asarray(y)
-        y2 = y * y / dt
-        return cls(y2, y2)
-
-    @classmethod
-    def pk(cls, y: np.ndarray, pk: np.ndarray, dt: float = 1.0) -> "GarchEstimator":
-        y = np.asarray(y)
-        y2 = y * y / dt
-        return cls(y2, np.asarray(pk) / dt)
-
-    def filter(self, p: np.ndarray) -> np.ndarray:
-        w, a, b = p
-        sig2 = np.zeros(self.n)
-        for i in range(self.n):
-            if i == 0:
-                sig2[0] = w / (1 - a - b)
-            else:
-                sig2[i] = w + a * self.p[i - 1] + b * sig2[i - 1]
-        return sig2
-
-    def log_like(self, p: np.ndarray) -> float:
-        """Log-likelihood of the GARCH model
-
-        The sign is flipped because we want to maximize not minimize
-        """
-        sig2 = self.filter(p)
-        return np.sum(np.log(sig2) + self.y2 / sig2)
-
-    def fit(self, **options: Any) -> dict:
-        p = np.array([0.01, 0.05, 0.94])
-        r = minimize(self.log_like, p, bounds=((0.001, None),) * 3, options=options)
-        if r.success:
-            return dict(params=r.x, aic=akaike_information(r.x, -r.fun))
-        raise RuntimeError(r)
diff --git a/quantflow_tests/test_utils.py b/quantflow_tests/test_utils.py
new file mode 100644
index 00000000..88c734f9
--- /dev/null
+++ b/quantflow_tests/test_utils.py
@@ -0,0 +1,40 @@
+import numpy as np
+
+from quantflow.utils.numbers import round_to_step, to_decimal
+from quantflow.utils.paths import Paths
+
+
+def test_round_to_step():
+    assert str(round_to_step(1.234, 0.1)) == "1.2"
+    assert str(round_to_step(1.234, 0.01)) == "1.23"
+    assert str(round_to_step(1.236, 0.01)) == "1.24"
+    assert str(round_to_step(1.1, 0.01)) == "1.10"
+    assert str(round_to_step(1.1, 0.001)) == "1.100"
+    assert str(round_to_step(2, 0.001)) == "2.000"
+    assert str(round_to_step(to_decimal("2.00000000000"), 0.001)) == "2.000"
+
+
+def test_normal_draws() -> None:
+    paths = Paths.normal_draws(100, 1, 1000)
+    assert paths.samples == 100
+    assert paths.time_steps == 1000
+    m = paths.mean()
+    np.testing.assert_array_almost_equal(m, 0)
+    paths = Paths.normal_draws(100, 1, 1000, antithetic_variates=False)
+    assert np.abs(paths.mean().mean()) > np.abs(m.mean())
+
+
+def test_normal_draws1() -> None:
+    paths = Paths.normal_draws(1, 1, 1000)
+    assert paths.samples == 1
+    assert paths.time_steps == 1000
+    paths = Paths.normal_draws(1, 1, 1000, antithetic_variates=False)
+    assert paths.samples == 1
+    assert paths.time_steps == 1000
+
+
+def test_path_stats() -> None:
+    paths = Paths.normal_draws(paths=2, time_steps=1000)
+    assert paths.paths_mean().shape == (2,)
+    assert paths.paths_std(scaled=True).shape == (2,)
+    assert paths.paths_var(scaled=False).shape == (2,)
diff --git a/quantflow_tests/test_weiner.py b/quantflow_tests/test_weiner.py
index 23ad9ed8..e789c34c 100644
--- a/quantflow_tests/test_weiner.py
+++ b/quantflow_tests/test_weiner.py
@@ -2,7 +2,6 @@
 import pytest
 
 from quantflow.sp.weiner import WeinerProcess
-from quantflow.utils.paths import Paths
 from quantflow_tests.utils import characteristic_tests
 
 
@@ -33,25 +32,6 @@ def test_sampling(weiner: WeinerProcess) -> None:
     assert std[0] == 0
 
 
-def test_normal_draws() -> None:
-    paths = Paths.normal_draws(100, 1, 1000)
-    assert paths.samples == 100
-    assert paths.time_steps == 1000
-    m = paths.mean()
-    np.testing.assert_array_almost_equal(m, 0)
-    paths = Paths.normal_draws(100, 1, 1000, antithetic_variates=False)
-    assert np.abs(paths.mean().mean()) > np.abs(m.mean())
-
-
-def test_normal_draws1() -> None:
-    paths = Paths.normal_draws(1, 1, 1000)
-    assert paths.samples == 1
-    assert paths.time_steps == 1000
-    paths = Paths.normal_draws(1, 1, 1000, antithetic_variates=False)
-    assert paths.samples == 1
-    assert paths.time_steps == 1000
-
-
 def test_support(weiner: WeinerProcess) -> None:
     m = weiner.marginal(0.01)
     pdf = m.pdf_from_characteristic(32)