Vasicek with jumps

rl.SwitchingVasicekJumps(chain, r0, a, b, sigma, jumpIntensity, jumpMean): the mean level, the volatility and the jump intensity may switch. QuantLib has no jump short-rate model; the frozen limit is checked against the affine closed form.

The model

\[dr_t = a\,(b_{y_t} - r_t)\,dt + \sigma_{y_t}\,dW_t + dJ_t,\]

where, while \(y_t = i\), \(J\) jumps at rate \(\ell_i\) by exponentially distributed amounts with mean \(m\).

Reduction

The pricing equation is the Vasicek one with the extra term \(\ell_i\,\mathbb{E}[u_i(t, r + J) - u_i(t, r)]\). With \(u_i = e^{-B(t)r}a_i(t)\) a jump of size \(J\) multiplies \(e^{-Br}\) by \(e^{-BJ}\), and for exponential jumps \(\mathbb{E}[e^{-BJ}] = 1/(1 + mB)\). So the jumps change only the forcing:

\[P_i(0, T) = e^{-B(T)\,r_0}\,a_i(T), \qquad g_i(t) = -a\,b_i\,B + \tfrac12\sigma_i^2\,B^2 + \ell_i\Big(\frac{1}{1 + mB} - 1\Big).\]

Closed form

Two integrals beyond the \(I_n\) of the Vasicek page are needed,

\[J_1 = \int_0^T\frac{dt}{1 + mB} = \frac{a}{a + m}\Big(T + \frac1a\log\Big(1 + \frac ma(1 - E)\Big)\Big), \qquad J_2 = \int_0^T\frac{dt}{(1 + mB)^2} = \frac{a}{a + m}\Big(J_1 - \frac1a\Big(\frac{1}{1 + \frac ma(1 - E)} - 1\Big)\Big),\]

with \(E = e^{-aT}\). Then, with \(\bar\ell, \tilde\ell = \tfrac12(\ell_1 \pm \ell_2)\),

\[\int_0^T\bar g = -a\bar b\,I_1 + \tfrac12\bar s\,I_2 + \bar\ell\,(J_1 - T),\]

\[\int_0^T\tilde g^{\,2} = a^2\tilde b^2 I_2 - a\tilde b\,\tilde s\,I_3 + \tfrac14\tilde s^2 I_4 + \tilde\ell^{\,2}(J_2 - 2J_1 + T) - 2a\tilde b\,\tilde\ell\Big(\frac{T - J_1}{m} - I_1\Big) + \tilde s\,\tilde\ell\Big(\frac1m\Big(I_1 - \frac{T - J_1}{m}\Big) - I_2\Big),\]

\[\tilde g(T) = -a\tilde b\,B + \tfrac12\tilde s\,B^2 + \tilde\ell\Big(\frac{1}{1 + mB} - 1\Big), \qquad \tilde g(0) = 0, \qquad \tilde g'(T) = \Big(-a\tilde b + \tilde s\,B - \frac{\tilde\ell\,m}{(1 + mB)^2}\Big)E.\]

Python

r0, a, b, sigma, ell, m, lam, T = 0.0, 2.0, [0.05, 0.02], [0.02, 0.01], [3.0, 0.2], 0.03, 10.0, 3.0

E, B, I1, I2, I3, I4 = powers_of_B(a, T)
J1 = a / (a + m) * (T + np.log(1 + m / a * (1 - E)) / a)
J2 = a / (a + m) * (J1 - (1 / (1 + m / a * (1 - E)) - 1) / a)
bbar, bt = half(b)
sbar, st = half([s * s for s in sigma])
lbar, lt = half(ell)
closed = np.exp(-B * r0) * second_order(
    int_gbar=-a * bbar * I1 + sbar * I2 / 2 + lbar * (J1 - T),
    int_gt2=(a * a * bt * bt * I2 - a * bt * st * I3 + st * st * I4 / 4 + lt * lt * (J2 - 2 * J1 + T)
             - 2 * a * bt * lt * ((T - J1) / m - I1) + st * lt * ((I1 - (T - J1) / m) / m - I2)),
    gt_T=-a * bt * B + st * B * B / 2 + lt * (1 / (1 + m * B) - 1),
    gt_0=0.0,
    dgt_T=(-a * bt + st * B - lt * m / (1 + m * B) ** 2) * E,
    lam=lam)

model = rl.SwitchingVasicekJumps(rl.RegimeChain.twoState(lam, lam), r0, a, b, sigma, ell, m)
print(f"closed form, second order  {closed:.8f}")
print(f"library, numerical         {library_bond(model, T):.8f}")

closed form, second order  0.86213677
library, numerical         0.86213670

regimelib.symbolic.VasicekJumpsBondFirstOrder gives the first-order formula for any finite chain as a sympy expression (Symbolic).