Other instruments
How each remaining instrument is priced from the reduced system \(a' = (Q + \operatorname{diag} g)\,a\): what is solved, with which terminal vector, and what is then integrated.
Bonds and claims on the regime
A bond, or a survival probability when the state is a default intensity, is \(P_i(0, T) = e^{-B(T)\cdot x_0}a_i(T)\) with \(a(0) = \mathbf 1\). The terminal vector says what is paid in each regime at \(T\). Starting from the unit vector \(e_j\) gives the value of one unit paid only if the regime at \(T\) is \(j\), and these add up to the bond.
lam, T = 5.0, 5.0
model = rl.SwitchingVasicek(rl.RegimeChain.twoState(lam, lam), 0.03, 0.5, [0.06, 0.02], [0.015, 0.008])
engine = rl.NumericalSwitchingEngine(model, regime=0)
parts = []
for j in (0, 1):
claim = rl.ZeroCouponBond(T, regimeAtMaturity=j)
claim.setPricingEngine(engine)
parts.append(claim.NPV())
print(f"paid only if the regime at T is {j} {parts[-1]:.8f}")
print(f"sum {sum(parts):.8f}")
print(f"bond {library_bond(model, T):.8f}")
paid only if the regime at T is 0 0.41647096
paid only if the regime at T is 1 0.41655024
sum 0.83302120
bond 0.83302120
European and digital options on equity
With \(\phi_i(u)\) the characteristic function of the martingale log return, a call is Lewis’s formula (the
lewis_call helper), and a digital call paying one if \(S_T > K\) is the Gil–Pelaez inversion of the same
function,
Sensitivities to the spot follow by differentiating under the integral. The derivative in \(T\) needs \(\partial_T\phi\), which the reduced system supplies without a further solve: \(\partial_Ta(T) = (Q + \operatorname{diag} g(T))\,a(T)\).
S0, r, q, sigma, lam, T, K = 100.0, 0.03, 0.0, [0.30, 0.15], 25.0, 1.0, 105.0
F = S0 * np.exp((r - q) * T)
def phi(u):
g = [-0.5 * s * s * (u * u + 1j * u) for s in sigma]
return exact_two_state(g[0], g[1], lam, T)
x, w = np.polynomial.legendre.leggauss(400)
u, w = (x + 1) * 30.0, w * 30.0 # nodes on (0, 60)
closed = np.exp(-r * T) * (0.5 + np.sum(w * (np.exp(1j * u * np.log(F / K)) * phi(u) / (1j * u)).real) / np.pi)
model = rl.SwitchingBlackScholesProcess(rl.RegimeChain.twoState(lam, lam), S0, r, q, sigma)
digital = rl.VanillaOption(("cash", "call", K, 1.0), maturity=T)
digital.setPricingEngine(rl.NumericalSwitchingEngine(model, regime=0))
print(f"digital call, closed form {closed:.6f} library {digital.NPV():.6f}")
digital call, closed form 0.408966 library 0.408966
Credit default swaps
When the state is a default intensity, \(Q(t) = e^{-B(t)x_0}a_i(t)\) is the survival probability. A swap paying the spread \(s\) on dates \(t_1 < \dots < t_n\), with accrued premium paid at default and recovery \(R\), has
with \(\tau_k = t_k - t_{k-1}\), \(t_k^m\) the midpoint of the period and \(D\) the discount curve. The fair spread is the ratio of the protection leg to the premium leg per unit spread. Only survival probabilities are needed, so the closed form is the Vasicek bond evaluated at each date.
x0, a, b, sigma, lam = 0.02, 0.5, [0.06, 0.01], [0.01, 0.004], 5.0
recovery, rate, times = 0.4, 0.03, [0.25 * k for k in range(1, 21)]
bbar, bt = half(b)
sbar, st = half([s * s for s in sigma])
def survival(t):
if t == 0:
return 1.0
E, B, I1, I2, I3, I4 = powers_of_B(a, t)
return np.exp(-B * x0) * second_order(
-a * bbar * I1 + sbar * I2 / 2, a * a * bt * bt * I2 - a * bt * st * I3 + st * st * I4 / 4,
-a * bt * B + st * B * B / 2, 0.0, (-a * bt + st * B) * E, lam)
premium = protection = 0.0
for t0, t1 in zip([0.0] + times[:-1], times):
tau, mid, dQ = t1 - t0, (t0 + t1) / 2, survival(t0) - survival(t1)
premium += tau * np.exp(-rate * t1) * survival(t1) + 0.5 * tau * np.exp(-rate * mid) * dQ
protection += (1 - recovery) * np.exp(-rate * mid) * dQ
model = rl.SwitchingVasicek(rl.RegimeChain.twoState(lam, lam), x0, a, b, sigma)
cds = rl.CreditDefaultSwap("buyer", 0.01, times, recovery, discount=rate)
cds.setPricingEngine(rl.NumericalSwitchingEngine(model, regime=0))
print(f"fair spread, closed form second order {1e4 * protection / premium:.4f} bp")
print(f"fair spread, library {1e4 * cds.fairSpread():.4f} bp")
fair spread, closed form second order 178.1448 bp
fair spread, library 178.1437 bp
Several names on one chain have joint survival equal to the bond of the summed forcing, so a first-to-default swap is a single-name swap on the sum (Intensity basket).
Geometric Asian options
Under Black–Scholes with a switching volatility the continuous geometric average is \(G = S_0e^Y\) with
Each increment of \(\log S\) enters \(Y\) scaled by the fraction of the life that remains, which in time to maturity \(\tau\) is \(\tau/T\). So \(\mathbb{E}[e^{iuY} \mid y_0 = i] = a_i(T)\) with
a polynomial in \(\tau\), so every integral is elementary:
The first integral is the Kemna–Vorst result: the average of a lognormal path has half the drift and a third of the variance. The option is Lewis’s formula on \(G\), whose forward is \(F_G = S_0\,\mathbb{E}[e^Y]\).
S0, r, q, sigma, lam, T = 100.0, 0.03, 0.0, [0.30, 0.15], 50.0, 1.0
sbar, st = half([s * s for s in sigma])
def phi_Y(u):
return second_order(
int_gbar=0.5j * u * (r - q - sbar / 2) * T - u * u * sbar * T / 6,
int_gt2=st * st * T / 4 * (-u * u / 3 + 1j * u ** 3 / 2 + u ** 4 / 5),
gt_T=-st / 2 * (1j * u + u * u),
gt_0=0.0,
dgt_T=-st / (2 * T) * (1j * u + 2 * u * u),
lam=lam)
def asian_call(K, U=60.0, n=400):
FG = S0 * phi_Y(-1j).real
k, lf = np.log(FG / K), np.log(FG / S0)
x, w = np.polynomial.legendre.leggauss(n)
u, w = (x + 1) * U / 2, w * U / 2
z = u - 0.5j
integral = np.sum(w * (np.exp(1j * u * k) * phi_Y(z) * np.exp(-1j * z * lf)).real / (u * u + 0.25))
return np.exp(-r * T) * (FG - np.sqrt(FG * K) / np.pi * integral)
model = rl.SwitchingBlackScholesProcess(rl.RegimeChain.twoState(lam, lam), S0, r, q, sigma)
for K in (95.0, 105.0):
asian = rl.ContinuousGeometricAsianOption(("call", K), maturity=T)
asian.setPricingEngine(rl.NumericalSwitchingEngine(model, regime=0))
print(f"K = {K:5.0f} closed form, second order {asian_call(K):.6f} library {asian.NPV():.6f}")
K = 95 closed form, second order 8.672621 library 8.672626
K = 105 closed form, second order 3.785976 library 3.785979
Barrier, American and Bermudan instruments
These depend on the path of the state and not only on its terminal value, so the state does not factor out and there is no reduced system. They are priced on the coupled equations themselves, one unknown function per regime,
on a grid in the state with one block per regime. A knock-out barrier is the condition \(u_i = \text{rebate}\) at the barrier in every regime, and a knock-in is the vanilla less the knock-out. An American option is the variational inequality
with an exercise boundary that differs by regime. A Bermudan swaption projects onto the exercise value at each date, and that value is the coupon bond in each regime from the reduced system.
lam = 5.0
model = rl.SwitchingBlackScholesProcess(rl.RegimeChain.twoState(lam, lam), 100.0, 0.05, 0.0, [0.30, 0.15])
grid = rl.SwitchingFDEngine(model, regime=0, n=801, steps=400)
european = rl.VanillaOption(("put", 100.0), maturity=1.0)
american = rl.VanillaOption(("put", 100.0), exercise="american", maturity=1.0)
knock_out = rl.BarrierOption("downout", 80.0, 0.0, ("put", 100.0), maturity=1.0)
knock_in = rl.BarrierOption("downin", 80.0, 0.0, ("put", 100.0), maturity=1.0)
for name, instrument in (("European put", european), ("American put", american),
("down-and-out put", knock_out), ("down-and-in put", knock_in)):
instrument.setPricingEngine(grid)
print(f"{name:18s} {instrument.NPV():.4f}")
print(f"{'out + in':18s} {knock_out.NPV() + knock_in.NPV():.4f}")
European put 7.1739
American put 7.7214
down-and-out put 1.1979
down-and-in put 5.9760
out + in 7.1739
Sensitivities to the parameters
The closed forms are explicit in the parameters, so their derivatives are formulas. To first order the log price is \(\int\bar g + \tfrac12\varepsilon\int\tilde g^{\,2} + \log(1 \pm \tfrac12\varepsilon\tilde g(T))\). If a parameter \(p_i\) of regime \(i\) enters the forcing linearly, \(g_i = p_i\,f(t) + \dots\), then
The first term is the sensitivity of the averaged model, shared equally between the regimes; the second moves it towards the regime whose parameter matters more; the third is the memory of the start. For the Vasicek mean level \(f = -aB\).
x0, a, b, sigma, lam, T = 0.02, 0.5, [0.06, 0.01], [0.01, 0.004], 5.0, 5.0
E, B, I1, I2, I3, I4 = powers_of_B(a, T)
_, bt = half(b)
_, st = half([s * s for s in sigma])
eps = 1 / lam
cross = -a * bt * I2 + st * I3 / 2 # int g~ B
formula = {"b1": -a / 2 * I1 - eps * a / 2 * cross - eps * a / 4 * B,
"b2": -a / 2 * I1 + eps * a / 2 * cross + eps * a / 4 * B,
"lambda": -eps ** 2 * (0.5 * (a * a * bt * bt * I2 - a * bt * st * I3 + st * st * I4 / 4)
+ 0.5 * (-a * bt * B + st * B * B / 2))}
def log_bond(b_, lam_):
return np.log(library_bond(rl.SwitchingVasicek(rl.RegimeChain.twoState(lam_, lam_), x0, a, b_, sigma), T))
d = 1e-5
numeric = {"b1": (log_bond([b[0] + d, b[1]], lam) - log_bond([b[0] - d, b[1]], lam)) / (2 * d),
"b2": (log_bond([b[0], b[1] + d], lam) - log_bond([b[0], b[1] - d], lam)) / (2 * d),
"lambda": (log_bond(b, lam + 1e-3) - log_bond(b, lam - 1e-3)) / 2e-3}
for name in formula:
print(f"d log P / d {name:7s} formula {formula[name]:+.6f} finite difference of the library {numeric[name]:+.6f}")
d log P / d b1 formula -1.622192 finite difference of the library -1.622189
d log P / d b2 formula -1.541978 finite difference of the library -1.541981
d log P / d lambda formula +0.000429 finite difference of the library +0.000426
For a general chain the first-order price depends on the chain through the stationary distribution \(\pi\) and the group inverse \(Q^\#\), whose derivatives are exact:
parameterGreek on the first-order symbolic classes implements this (Symbolic).