Variance gamma

rl.SwitchingVarianceGammaProcess(chain, S0, r, q, sigma, nu, theta) is QuantLib’s VarianceGammaProcess with all three parameters allowed to switch.

The model

While \(y_t = i\) the log price is a Brownian motion with drift \(\theta_i\) and volatility \(\sigma_i\) run on a gamma clock with variance rate \(\nu_i\):

\[d\log S_t = (r - q + \omega_i)\,dt + dX^{(i)}_t, \qquad X^{(i)}_t = \theta_i\,G^{(i)}_t + \sigma_i\,W_{G^{(i)}_t}, \qquad \omega_i = \frac{1}{\nu_i}\log\big(1 - \theta_i\nu_i - \tfrac12\sigma_i^2\nu_i\big).\]

Reduction

Given the clock increment \(dG\), the increment of \(X\) is normal with mean \(\theta_i\,dG\) and variance \(\sigma_i^2\,dG\), and the gamma increment has \(\mathbb{E}[e^{-z\,dG}] = (1 + \nu_iz)^{-dt/\nu_i}\). Putting \(z = -iu\theta_i + \tfrac12\sigma_i^2u^2\) gives the Lévy exponent, and for the martingale log return

\[\phi_i(u) = a_i(T), \qquad g_i = \psi_i(u) + iu\,\omega_i, \qquad \psi_i(u) = -\frac{1}{\nu_i}\log\big(1 - iu\,\theta_i\nu_i + \tfrac12\sigma_i^2\nu_i\,u^2\big).\]

All three parameters sit inside \(g_i\) and there is no state, hence no Riccati equation for them to disturb. That is why all three may switch.

Closed form

The forcing is constant, so the two-regime characteristic function is exact, as for Black–Scholes, with \(\bar g, \tilde g = \tfrac12(g_1 \pm g_2)\). The averaged model, \(e^{\bar gT}\), is the sum of two independent variance gamma processes each run at half speed; it is not itself variance gamma unless the regimes agree. The characteristic function decays like a power of \(u\), so the Fourier integral is carried further than for a diffusion.

Python

S0, r, q, sigma, nu, theta, lam, T = 100.0, 0.03, 0.0, [0.25, 0.12], [0.5, 0.2], [-0.25, -0.10], 25.0, 1.0

def phi(u):
    g = []
    for s, n, th in zip(sigma, nu, theta):
        omega = np.log(1 - th * n - 0.5 * s * s * n) / n
        g.append(-np.log(1 - 1j * u * th * n + 0.5 * s * s * n * u * u) / n + 1j * u * omega)
    return exact_two_state(g[0], g[1], lam, T)

model = rl.SwitchingVarianceGammaProcess(rl.RegimeChain.twoState(lam, lam), S0, r, q, sigma, nu, theta)
for K in (90.0, 110.0):
    print(f"K = {K:5.0f}   closed form, exact  {lewis_call(phi, S0, K, r, q, T, U=600.0, n=6000):.6f}"
          f"   library  {library_call(model, K, T):.6f}")

K =    90   closed form, exact  16.402558   library  16.402558
K =   110   closed form, exact  5.107950   library  5.107950