BNS Model

The Barndorff-Nielson–Shephard (BNS) model is a stochastic volatility model where the variance process \(\nu_t\), or better, the activity rate process, follows a non-gaussian OU process. The leverage effect can be accommodated by correlating the Brownian motion \(w_t\) and the BDLP \(z_t\) as the following equations illustrate:

(63)\[\begin{align} y_t &= w_{\tau_t} + \rho z_{\kappa t} \\ d \nu_t &= -\kappa \nu_t dt + d z_{\kappa t} \\ \tau_t &= \int_0^t \nu_s ds \end{align}\]

This means that the characteristic function of \(y_t\) can be represented as

(64)\[\begin{align} \Phi_{y_t, u} & = {\mathbb E}\left[\exp{\left(i u w_{\tau_t} + i u \rho z_{\kappa t}\right)}\right] \\ &= {\mathbb E}\left[\exp{\left(-\tau_t \phi_{w, u} + i u \rho z_{\kappa t}\right)}\right] \end{align}\]

\(\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 process section.

from quantflow.sp.bns import BNS

pr = BNS.create(vol=0.5, decay=10, kappa=10, rho=-1)
pr

BNS(variance_process=GammaOU(rate=0.25, kappa=10.0, bdlp=CompoundPoissonProcess[Exponential](intensity=2.5, jumps=Exponential(decay=10.0))), rho=-1.0)

from quantflow.utils import plot
m = pr.marginal(2)
plot.plot_characteristic(m, max_frequency=10)

Marginal Distribution

m.mean(), m.std()

(1.13220528467725e-08, 0.524669184061491)

plot.plot_marginal_pdf(m, 128, normal=True, analytical=False)

Appendix

Carr at al [CPM03] show that the join characteristic function of \(\tau_t\) and \(z_{\kappa t}\) has a closed formula, and this is our derivation

(65)\[\begin{align} \zeta_{a, b} &= \ln {\mathbb E} \left[\exp{\left(i a \tau_t + i b z_{\kappa t}\right)}\right] \\ \zeta_{a, b} &= i c \nu_0 - \int_b^{b+c} \frac{\phi_{z_1, s}}{a+\kappa b - \kappa s} ds = i c \nu_0 + \lambda \left(I_{b+c} - I_{b}\right) \\ c &= a \frac{1 - e^{-\kappa t}}{\kappa} \end{align}\]

Noting that (see non-gaussian OU process)

(66)\[\begin{align} i a \tau_t + i b z_{\kappa t} &= i a \epsilon_t \nu_0 + \int_0^t \left(i a \epsilon_{t-s} + i b\right) d z_{\kappa s} \\ &= i a \epsilon_t \nu_0 + \int_0^{\kappa t} \left(i a \epsilon_{t-s/\kappa} + i b\right) d z_s \\ \epsilon_t &= \frac{1 - e^{-\kappa t}}{\kappa} \end{align}\]

we obtain

(67)\[\begin{align} \zeta_{a, b} &= i a \epsilon_t \nu_0 + \ln {\mathbb E} \left[\exp{\left(\int_0^{\kappa t} \left(i a \epsilon_{t-s/\kappa} + i b\right) d z_s\right)}\right] \\ &= i a \epsilon_t \nu_0 - \int_0^{\kappa t} \phi_z\left(a \epsilon_{t-s/\kappa} + b\right) d s \\ &= i a \epsilon_t \nu_0 - \int_L^U \frac{\phi_{z,s}}{a + \kappa b - \kappa s} d s \end{align}\]

Here we use sympy to derive the integral in the characteristic function.

import sympy as sym

k =  sym.Symbol("k")
iβ = sym.Symbol("iβ")
γ = sym.Symbol("γ")
s = sym.Symbol("s")
ϕ = s/(s+iβ)/(γ-k*s)
ϕ

\[\displaystyle \frac{s}{\left(iβ + s\right) \left(- k s + γ\right)}\]

r = sym.integrate(ϕ, s)
sym.simplify(r)

\[\displaystyle \frac{- iβ k \log{\left(iβ + s \right)} - γ \log{\left(\frac{k s - γ}{k} \right)}}{k \left(iβ k + γ\right)}\]

import numpy as np
f = lambda x: x*np.log(x)
f(0.001)

-0.006907755278982137