Models: mathematics and Python

One page per model. Each page states the model, reduces its pricing problem to a linear system, gives the closed form for a two-regime chain, and then evaluates that closed form in a few lines of Python next to the library’s price for the same inputs.

Short rates and credit

Model	Class	Switches	Forcing \(g_i\)	Closed form
Vasicek	SwitchingVasicek	\(b,\ \sigma\)	\(-a b_i B + \tfrac12\sigma_i^2 B^2\)	second order
Cox–Ingersoll–Ross	SwitchingCoxIngersollRoss	\(\theta\)	\(-k\theta_i B\)	second order
Hull–White	SwitchingHullWhite	\(\sigma\)	\(\tfrac12\sigma_i^2 B^2\)	second order
G2++	SwitchingG2	\(\sigma,\ \eta,\ \rho\)	\(\tfrac12\sigma_i^2B_a^2 + \tfrac12\eta_i^2B_b^2 + \rho_i\sigma_i\eta_iB_aB_b\)	second order
Intensity basket	SwitchingIntensityBasket	each name’s	sum of the names’ forcings	first order, for the dependence

Equity diffusions

Model	Class	Switches	Forcing \(g_i\)	Closed form
Black–Scholes	SwitchingBlackScholesProcess	\(\sigma\)	\(-\tfrac12\sigma_i^2(u^2 + iu)\)	exact
Heston	SwitchingHestonModel	\(\theta\)	\(\kappa\theta_i D(t)\)	second order
CEV	SwitchingCEVProcess	\(\sigma\)	\(-z\sigma_i^2h(\tau)\)	exact mixture; second order
Heston with a switching volatility of variance	SwitchingHestonVolOfVol	\(\xi\)	none: the switched operators do not commute	first-order rule

Jump models

Model	Class	Switches	Forcing \(g_i\)	Closed form
Vasicek with jumps	SwitchingVasicekJumps	\(b,\ \sigma,\ \ell\)	\(-a b_i B + \tfrac12\sigma_i^2 B^2 + \ell_i\big((1 + mB)^{-1} - 1\big)\)	second order
Merton jump diffusion	SwitchingMerton76Process	\(\sigma,\ \ell\)	\(-\tfrac12\sigma_i^2(u^2 + iu) + \ell_i(\phi_J - 1 - iu\bar k)\)	exact
Bates	SwitchingBatesModel	\(\theta,\ \ell\)	\(\kappa\theta_i D(t) + \ell_i(\phi_J - 1 - iu\bar k)\)	second order
Variance gamma	SwitchingVarianceGammaProcess	\(\sigma,\ \nu,\ \theta\)	\(\psi_i(u) + iu\,\omega_i\)	exact

Hybrid

Model	Class	Switches	Forcing \(g_i\)	Closed form
Equity with stochastic rates	SwitchingEquityRates	both parts’	rates forcing at weight \(1 - iw\) plus equity forcing at \(w\)	second order

Instruments

Page	Covers
Options on bonds, swaptions and caps	options on zero-coupon and coupon bonds, swaptions, caps: the regime at expiry, Gil–Pelaez inversion, Jamshidian’s decomposition by regime, the first-order closed form
Other instruments	claims on the regime, digital options, credit default swaps, geometric Asian options, barrier and American options, sensitivities to the parameters
Observed and inferred regimes	a belief as the starting regime; observed against inferred regimes for options on bonds and swaps and for early exercise; the belief as a state variable

The common formula

A price under a switching model is one function per regime, \(u_i(t, x)\), and the pricing equations are coupled through the generator \(Q\) of the chain. For every model in the tables except Heston with a switching volatility of variance, the state factors out and what remains is a linear system in time alone,

\[a'(t) = \big(Q + \operatorname{diag} g(t)\big)\,a(t), \qquad a(0) = \mathbf 1,\]

with one forcing \(g_i(t)\) per regime. The model pages differ only in \(g_i\) and in the prefactor that multiplies \(a_i\).

The closed forms are written for two regimes switching at rate \(\lambda\) in each direction, rl.RegimeChain.twoState(lam, lam), with

\[\varepsilon = \frac1\lambda, \qquad \bar g = \tfrac12(g_1 + g_2), \qquad \tilde g = \tfrac12(g_1 - g_2).\]

The average \(\bar g\) is the forcing of the averaged model, and the half-difference \(\tilde g\) measures how much the regimes disagree. To second order in \(\varepsilon\), up to terms of size \(e^{-2\lambda T}\),

\[a_{1,2}(T) = \exp\Big(\int_0^T \bar g + \frac{\varepsilon}{2}\int_0^T \tilde g^{\,2} - \frac{\varepsilon^2}{8}\big(\tilde g(T)^2 + \tilde g(0)^2\big)\Big) \Big(1 \pm \frac{\varepsilon}{2}\,\tilde g(T) \mp \frac{\varepsilon^2}{4}\,\tilde g'(T)\Big) + O(\varepsilon^3),\]

with the upper sign for a start in the first regime (regime=0). The three parts have fixed meanings:

• \(e^{\int\bar g}\) is the averaged model;

• \(\tfrac12\varepsilon\int\tilde g^{\,2}\) is the Green–Kubo term, half the variance of the fluctuating integral of \(g\), the same from both regimes;

• the bracket is the memory of the starting regime.

When \(g\) does not depend on time the system is solved exactly:

\[a_{1,2}(T) = e^{(\bar g - \lambda)T}\Big(\cosh sT + \frac{\lambda \pm \tilde g}{s}\,\sinh sT\Big), \qquad s = \sqrt{\lambda^2 + \tilde g^2}.\]

So each page needs five things from its model: \(\int\bar g\), \(\int\tilde g^{\,2}\), \(\tilde g(T)\), \(\tilde g(0)\) and \(\tilde g'(T)\).

The helpers

The Python on every page uses these functions. They are the two formulas above, the integrals of \(B(t) = (1 - e^{-at})/a\), Lewis’s formula, and two wrappers around the library.

import numpy as np
import regimelib as rl


def second_order(int_gbar, int_gt2, gt_T, gt_0, dgt_T, lam, sign=+1):
    """a_i(T) for a symmetric two-state chain through eps^2 = 1 / lam^2. sign=+1 is a start in regime 0."""
    eps = 1.0 / lam
    return (np.exp(int_gbar + eps / 2 * int_gt2 - eps ** 2 / 8 * (gt_T ** 2 + gt_0 ** 2))
            * (1 + sign * eps / 2 * gt_T - sign * eps ** 2 / 4 * dgt_T))


def exact_two_state(g1, g2, lam, T, sign=+1):
    """a_i(T) when the forcing is constant in time: exact, no expansion."""
    gbar, gt = (g1 + g2) / 2, (g1 - g2) / 2
    s = np.sqrt(lam ** 2 + gt ** 2)
    return np.exp((gbar - lam) * T) * (np.cosh(s * T) + (lam + sign * gt) / s * np.sinh(s * T))


def half(x):
    """Half-sum and half-difference across the two regimes."""
    return (x[0] + x[1]) / 2, (x[0] - x[1]) / 2


def powers_of_B(a, T):
    """E = e^{-aT}, B = (1 - E) / a and I_n = int_0^T B(t)^n dt for n = 1..4."""
    E = np.exp(-a * T)
    I1 = (T - (1 - E) / a) / a
    I2 = (T - 2 * (1 - E) / a + (1 - E ** 2) / (2 * a)) / a ** 2
    I3 = (T - 3 * (1 - E) / a + 3 * (1 - E ** 2) / (2 * a) - (1 - E ** 3) / (3 * a)) / a ** 3
    I4 = (T - 4 * (1 - E) / a + 3 * (1 - E ** 2) / a - 4 * (1 - E ** 3) / (3 * a) + (1 - E ** 4) / (4 * a)) / a ** 4
    return E, (1 - E) / a, I1, I2, I3, I4


def lewis_call(phi, S0, K, r, q, T, U=60.0, n=400):
    """European call from the characteristic function phi(u) of the martingale log return X, S_T = F e^X."""
    F = S0 * np.exp((r - q) * T)
    k = np.log(F / K)
    x, w = np.polynomial.legendre.leggauss(n)
    u, w = (x + 1) * U / 2, w * U / 2
    integral = np.sum(w * (np.exp(1j * u * k) * phi(u - 0.5j)).real / (u * u + 0.25))
    return np.exp(-r * T) * (F - np.sqrt(F * K) / np.pi * integral)


def library_bond(model, T, regime=0):
    bond = rl.ZeroCouponBond(T)
    bond.setPricingEngine(rl.NumericalSwitchingEngine(model, regime=regime))
    return bond.NPV()


def library_call(model, K, T, regime=0):
    option = rl.VanillaOption(("call", K), maturity=T)
    option.setPricingEngine(rl.NumericalSwitchingEngine(model, regime=regime))
    return option.NPV()

For a chain with more regimes or unequal rates the library’s engines apply as they stand: FastSwitchingEngine runs the same expansion to any order with the group inverse of \(Q\) in place of the division by \(2\lambda\), and NumericalSwitchingEngine solves the linear system.

Adding jumps

A regime switch is already a jump, in a parameter: the volatility, the level a rate reverts to, or the long-run variance moves at once to a new value, while the stock price or the rate itself stays continuous. Two further kinds of jump move the price or the rate itself, and both keep a closed form.

The price or the rate jumps at random times, more often in one regime than the other. Examples: a stock price that drops by a lognormal factor at Poisson times, twice a year in a stressed regime and once in five years in a calm one (Merton jump diffusion, and Bates when the variance is also stochastic); a short rate or a default intensity that jumps upward by an exponential amount, three times a year in one regime and rarely in the other (Vasicek with jumps); a stock price that moves only by jumps, with its skew and kurtosis set by the regime (Variance gamma). In each case the jump adds one term to the forcing, the regime’s jump rate times the transform of the jump size less one:

\[\text{rate jumps by } J:\quad \ell_i\,\big(\mathbb{E}[e^{-B(t)J}] - 1\big), \qquad\qquad \text{log price jumps by } J:\quad \ell_i\,\big(\mathbb{E}[e^{iuJ}] - 1 - iu\bar k\big).\]

Any jump-size law with a known transform will do. The jump rate may switch, and so may the size law, by giving each regime its own transform. For a stock the term is constant in time, so the two-regime characteristic function is exact. The one restriction is that a jump rate proportional to the rate or the variance itself enters the Riccati equation, and then only its constant part may switch.

The price or the rate jumps at the moment the regime changes. Examples: the stock falls five percent on entry to the crisis regime and recovers part of it on exit; the short rate steps down fifty basis points when the regime turns to easing. If a move from regime \(i\) to \(j\) shifts the log price or the rate by a random amount \(J_{ij}\), each off-diagonal switching rate is multiplied by the transform of that jump, \(\Phi_{ij} = \mathbb{E}[e^{iuJ_{ij}}]\) for a log price and \(\mathbb{E}[e^{-B(t)J_{ij}}]\) for a rate, and the system is \(a' = (Q\circ\Phi + \operatorname{diag} g)\,a\). It is still linear, and for a stock with two regimes it is still exact:

\[\phi_1(u) = e^{(\bar g - \lambda)T}\Big(\cosh sT + \frac{\tilde g + \lambda\Phi_{12}}{s}\,\sinh sT\Big), \qquad s = \sqrt{\tilde g^2 + \lambda^2\,\Phi_{12}\Phi_{21}},\]

with \(\phi_2\) given by exchanging the regimes. What changes is the expansion: these jumps arrive at the switching rate, so as the chain speeds up they pile up unless their size shrinks with the holding time, and the averaged model is no longer the model at the averaged parameters. The library implements the first kind and not this one.