From Lando (2004)* I am trying to replicate the following figure (Section 2.6 Default Barriers: The Black-Cox Setup):

The spreads are computed as follows:

$$s(T) = \frac{1}{T}\ln\frac{D}{B_0}-r$$

where:

• $$T$$ is the time to maturity of the bond,
• $$D$$ is the face value of the debt,
• $$B_0$$ is the bond price at time $$t=0$$,
• $$r$$ is the risk-free rate.

In the Black-Cox model, the bond price $$B_0$$ is given by:

$$B_0 = E^{\mathbb Q}\left[e^{-r\,T} \min\left(V_T, D\right)\mathbf 1_{\{\tau>T\}} \right] + E^{\mathbb Q}\left[e^{-r\,\tau} C_1(\tau) \mathbf 1_{\{\tau\leq T\}} \right]$$

where $$\tau$$ is the firs-hitting time of the barrier $$C_1(t) = e^{-\gamma\,(T-t)}C$$.

Analytical solution: $$B_0(V, T, D, C, \gamma, a) = BL(e^{-\gamma \,T} V, T, e^{-\gamma \,T} D, e^{-\gamma \,T} C) + B_b(V, T, C, \gamma, a)$$

with:

where $$H(V, T, L) = e^{-r\,T}N\left(\frac{\ln(V/L) + \left(r-\frac{1}{2}\sigma^2\right)T}{\sigma\sqrt{T}}\right)$$

and

where

$$b = \frac{\ln(C/V)-\gamma\,T}{\sigma}$$ $$\mu = \frac{r-a-\frac{1}{2}\sigma^2-\gamma}{\sigma}$$ $$\tilde \mu = \sqrt{\mu^2+2r}$$

My attempt, in Python (Jupyter Lab):

# Imports
import numpy as np
from scipy import stats

import plotly.graph_objects as go
from ipywidgets import interact

# Paramters
V = 120
D = 100
r = 0.05
σ = 0.2
T = 10

C = 90
γ = 0.02

# Time to maturity
dt = 1/1000
τ = np.arange(dt, T+dt, dt)

# Normal CDF
N = lambda x: stats.norm.cdf(x)

H = lambda V, T, L: np.exp(-r*T) * N( (np.log(V/L) + (r-0.5*σ**2)*T) / (σ*np.sqrt(T)) )

# Black-Scholes Call price
def C_BS(V, K, T):
d1 = (np.log(V/K) + (r + 0.5*σ**2)*T ) / ( σ*np.sqrt(T) )
d2 = d1 - σ*np.sqrt(T)
return V*N(d1) - np.exp(-r*T)*K*N(d2)

def BL(V, T, D, L):
return L * H(V, T, L) - L * (L/V)**(2*r/σ**2-1) * H(L**2/V, T, L) + \
C_BS(V, L, T) - (L/V)**(2*r/σ**2-1) * C_BS(L**2/V, L, T) - \
C_BS(V, D, T) + (L/V)**(2*r/σ**2-1) * C_BS(L**2/V, D, T)

def Bb(V, T, C, γ, a):
b = (np.log(C/V) - γ*T) / σ
μ = (r - a - 0.5*σ**2 - γ) / σ
m = np.sqrt(μ**2 + 2*r)
return C*np.exp(b*(μ-m)) * ( N((b-m*T)/np.sqrt(T)) + np.exp(2*m*b)*N((b+m*T)/np.sqrt(T)) )

def Black_Cox(V, T, C=90, a=0.001):
return BL(V*np.exp(-γ*T), T, D*np.exp(-γ*T), C*np.exp(-γ*T)) + Bb(V, T, C, γ, a)

s = lambda Model: 10000 * (np.log(D/Model) / τ - r)

# Coordinates for plotting
x = τ
y_Black_Cox = s(Black_Cox(V, τ, C=C, a=0.001))
y_Meton = s(V - C_BS(V, D, τ))


My Plotly plot looks like this:

fig = go.FigureWidget()
fig.update_layout(title='<b>Figure 2.9.</b> Black-Cox vs. Meton',
xaxis_title='Time to maturity', yaxis_title='Yield spread (bps)',
height=600, showlegend=False)
@interact(V=(50, 250, 1), C=(0, 300, 0.01), a=(1e-3, 1, 0.01))
def update(V=120, C=90, a=0.001):
with fig.batch_update():
fig.data[0].y = s(Black_Cox(V, τ, C=C, a=a))  # Black-Cox
fig.data[1].y = s(V - C_BS(V, D, τ))          # Merton
fig


While the spreads derived from the Merton model look the same, the ones from the Black-Cox do not match the ones in the figure above. What am I missing?

UPDATE:

I have changed the main function to correct the risk-free rate, so that the following relation holds:

def Black_Cox(V, T, C=90, a=0):
return np.exp(γ*T)*BL(V*np.exp(-γ*T), T, D*np.exp(-γ*T), C*np.exp(-γ*T)) + Bb(V, T, C, γ, a)


The plot looks better now, but still it doesn't resemble Fig. 2.9:

I have also implemented the result from the original paper**:

def Black_Cox1(V, C, T, a=0):
y = C*np.exp(-γ*T) / V
θ = (r - a - γ + 0.5*σ**2) / σ**2
δ = (r - a - γ - 0.5*σ**2)**2 + 2*σ**2*(r-γ)
ξ = np.sqrt(δ) / σ**2
η = np.sqrt(δ - 2*σ**2*a) / σ**2

z1 = (np.log(V/D) + (r-a-0.5*σ**2)*T) / (σ*np.sqrt(T))
z2 = (np.log(V/D) + 2*np.log(y) + (r-a-0.5*σ**2)*T) / (σ*np.sqrt(T))
z3 = (np.log(D/V) - (r-a+0.5*σ**2)*T) / (σ*np.sqrt(T))
z4 = (np.log(V/D) + 2*np.log(y) + (r-a+0.5*σ**2)*T) / (σ*np.sqrt(T))
z5 = (np.log(y) + ξ*σ**2*T) / (σ*np.sqrt(T))
z6 = (np.log(y) - ξ*σ**2*T) / (σ*np.sqrt(T))
z7 = (np.log(y) + η*σ**2*T) / (σ*np.sqrt(T))
z8 = (np.log(y) - η*σ**2*T) / (σ*np.sqrt(T))

return D*np.exp(-r*T)*(N(z1) - y**(2*θ-2)*N(z2)) + \
V*np.exp(-a*T)*(N(z3) + y**(2*θ)*N(z4) + y**(θ+ξ)*np.exp(a*T)*N(z5) + \
y**(θ-ξ)*np.exp(a*T)*N(z6) - y**(θ-η)*N(z7) - y**(θ-η)*N(z8))

y_Black_Cox_paper = s(Black_Cox1(V, C, τ))


And the plot looks like this:

Letting $$C \rightarrow 0$$ will make the blue and green lines identical to the one obtained from the Merton model. This makes sense: the barrier won't be hit and the bondholder, just like in the Merton model, will receive $$\min\left(V_T,D\right)$$.

I tend to believe that the plot was produced with $$\gamma = 0.002$$ instead of $$0.02$$.

* Lando, D. (2004). Credit Risk Modeling: Theory and Applications.

** Black, F., & Cox, J. C. (1976). Valuing corporate securities: Some effects of bond indenture provisions. The Journal of Finance, 31(2), 351-367.