xjtulyc Fixed Income Analytics
One-line summary: Price bonds, bootstrap yield curves, compute duration and convexity, fit Nelson-Siegel models, and analyze credit spreads using pure Python with numpy and scipy.
When to Use This Skill
- When pricing fixed-rate bonds, floaters, or zero-coupon bonds
- When computing modified duration, convexity, and DV01
- When bootstrapping spot rates from par yields or bond prices
- When fitting Nelson-Siegel or Svensson yield curve models
- When modeling credit spreads and CDS pricing
- When computing carry, roll-down, and term premium
Trigger keywords: bond pricing, yield to maturity, duration, convexity, DV01, PVBP, yield curve, Nelson-Siegel, Svensson, spot rate, forward rate, credit spread, CDS, OAS, zero curve, discount factor
Background & Key Concepts
Bond Pricing
A coupon bond with face value $F$, coupon $C$, yield $y$, and $T$ periods:
$$ P = \sum_{t=1}^{T} \frac{C}{(1+y)^t} + \frac{F}{(1+y)^T} $$
Modified Duration and Convexity
$$ \text{Duration} = \frac{1}{P}\sum_{t=1}^T \frac{t \cdot CF_t}{(1+y)^t} $$
$$ \text{Convexity} = \frac{1}{P}\sum_{t=1}^T \frac{t(t+1) \cdot CF_t}{(1+y)^{t+2}} $$
Price change: $\Delta P \approx -D^* \cdot P \cdot \Delta y + \frac{1}{2} C \cdot P \cdot (\Delta y)^2$
Nelson-Siegel Model
$$ y(m) = \beta_0 + \beta_1 \frac{1-e^{-m/\lambda}}{m/\lambda} + \beta_2 \left(\frac{1-e^{-m/\lambda}}{m/\lambda} - e^{-m/\lambda}\right) $$
where $m$ = maturity, $\beta_0$ = long-run level, $\beta_1$ = slope, $\beta_2$ = curvature.
Environment Setup
Install Dependencies
pip install numpy>=1.24 scipy>=1.11 pandas>=2.0 matplotlib>=3.7
Verify Installation
import numpy as np
from scipy.optimize import brentq
# Simple bond pricing test
def bond_price(face, coupon, y, n):
"""Annual coupon bond price given yield."""
t = np.arange(1, n+1)
pv = coupon / (1+y)**t
return pv.sum() + face / (1+y)**n
price = bond_price(1000, 50, 0.05, 10)
print(f"10Y 5% coupon bond at y=5%: ${price:.4f}")
# Expected: 1000.00 (par bond)
Core Workflow
Step 1: Bond Pricing and Duration/Convexity
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.optimize import brentq
# ------------------------------------------------------------------ #
# Bond analytics: price, YTM, duration, convexity, DV01
# ------------------------------------------------------------------ #
class Bond:
"""Fixed-rate bond with annual coupon payments."""
def __init__(self, face, coupon_rate, maturity, settlement_price=None, ytm=None):
"""
Parameters
----------
face : float — face (par) value
coupon_rate : float — annual coupon rate (e.g., 0.05 for 5%)
maturity : int — years to maturity
settlement_price : float — dirty price (optional)
ytm : float — yield to maturity (optional)
"""
self.face = face
self.coupon = face * coupon_rate
self.T = maturity
self.cashflows = np.array([self.coupon] * maturity)
self.cashflows[-1] += face # Final payment = coupon + face
self.times = np.arange(1, maturity + 1, dtype=float)
if ytm is not None:
self.ytm = ytm
self.price = self._price_from_ytm(ytm)
elif settlement_price is not None:
self.price = settlement_price
self.ytm = self._ytm_from_price(settlement_price)
else:
raise ValueError("Provide either ytm or settlement_price")
def _price_from_ytm(self, y):
"""Compute dirty price from YTM."""
discount = (1 + y) ** (-self.times)
return np.dot(self.cashflows, discount)
def _ytm_from_price(self, P):
"""Find YTM via root-finding."""
f = lambda y: self._price_from_ytm(y) - P
return brentq(f, -0.5, 10.0)
@property
def modified_duration(self):
"""Modified duration (years)."""
discount = (1 + self.ytm) ** (-self.times)
pv_cf = self.cashflows * discount
macaulay = np.dot(self.times, pv_cf) / self.price
return macaulay / (1 + self.ytm)
@property
def convexity(self):
"""Convexity (years²)."""
discount = (1 + self.ytm) ** (-self.times)
pv_cf = self.cashflows * discount
return np.dot(self.times * (self.times + 1), pv_cf) / (self.price * (1 + self.ytm)**2)
@property
def dv01(self):
"""Dollar value of 1 basis point (per face)."""
return self.modified_duration * self.price * 0.0001
def price_at_yield(self, y):
return self._price_from_ytm(y)
def summary(self):
return {
'Price': round(self.price, 4),
'YTM': round(self.ytm * 100, 4),
'Mod Duration': round(self.modified_duration, 4),
'Convexity': round(self.convexity, 4),
'DV01': round(self.dv01, 4),
}
# ---- Example bonds --------------------------------------------- #
bonds = {
'2Y 5% par bond': Bond(1000, 0.05, 2, ytm=0.05),
'10Y 5% par bond': Bond(1000, 0.05, 10, ytm=0.05),
'30Y 5% par bond': Bond(1000, 0.05, 30, ytm=0.05),
'10Y 3% discount': Bond(1000, 0.03, 10, ytm=0.05),
'10Y 0% zero': Bond(1000, 0.00, 10, ytm=0.05),
}
print(f"{'Bond':22s} {'Price':>8s} {'YTM%':>6s} {'Dur':>6s} {'Conv':>8s} {'DV01':>7s}")
print("-" * 70)
for name, bond in bonds.items():
s = bond.summary()
print(f"{name:22s} {s['Price']:8.2f} {s['YTM']:6.2f} {s['Mod Duration']:6.3f} "
f"{s['Convexity']:8.4f} {s['DV01']:7.4f}")
# ---- Price-yield relationship ---------------------------------- #
bond_10y = bonds['10Y 5% par bond']
yields = np.linspace(0.01, 0.15, 200)
prices = [bond_10y.price_at_yield(y) for y in yields]
# Duration approximation
dy = np.array(yields) - bond_10y.ytm
p_linear = bond_10y.price * (1 - bond_10y.modified_duration * dy)
p_convex = bond_10y.price * (1 - bond_10y.modified_duration * dy
+ 0.5 * bond_10y.convexity * dy**2)
fig, ax = plt.subplots(figsize=(9, 5))
ax.plot(yields * 100, prices, 'b-', linewidth=2.5, label='Exact price')
ax.plot(yields * 100, p_linear, 'r--', linewidth=1.5, label='Duration approx (linear)')
ax.plot(yields * 100, p_convex, 'g-.', linewidth=1.5, label='Duration + Convexity')
ax.axvline(bond_10y.ytm * 100, color='gray', linestyle=':', linewidth=1)
ax.axhline(bond_10y.price, color='gray', linestyle=':', linewidth=1)
ax.set_xlabel("Yield (%)"); ax.set_ylabel("Price ($)")
ax.set_title("Price-Yield Relationship — 10Y 5% Bond")
ax.legend(); ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("bond_price_yield.png", dpi=150)
plt.show()
Step 2: Yield Curve Bootstrapping
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.interpolate import CubicSpline
# ------------------------------------------------------------------ #
# Bootstrap zero/spot curve from par yields (Treasury benchmark)
# ------------------------------------------------------------------ #
# Par yield curve (observed market data)
maturities = np.array([0.5, 1, 2, 3, 5, 7, 10, 20, 30])
par_yields = np.array([5.20, 5.10, 4.90, 4.75, 4.55, 4.45, 4.30, 4.40, 4.35]) / 100
# ---- Bootstrap spot rates --------------------------------------- #
spot_rates = np.zeros(len(maturities))
for i, (T, y_par) in enumerate(zip(maturities, par_yields)):
if T <= 1:
# Short end: spot ≈ par (no intermediate coupons)
spot_rates[i] = y_par
else:
# Price = 1 (par bond); coupon = y_par per period (semi-annual)
# Sum of discounted intermediate coupons using already-bootstrapped spots
n_periods = int(T * 2) # Semi-annual
coupon = y_par / 2 # Semi-annual coupon rate
# Interpolate spots for intermediate maturities
t_interp = np.arange(0.5, T, 0.5) # Intermediate dates
if i > 0:
cs = CubicSpline(maturities[:i], spot_rates[:i])
spots_interp = cs(t_interp)
else:
spots_interp = np.array([y_par])
# Present value of intermediate coupons
pv_coupons = sum(coupon / (1 + s/2)**(2*t)
for t, s in zip(t_interp, spots_interp))
# Solve for final spot rate
# 1 = pv_coupons + (1 + coupon) / (1 + s_T/2)^(2T)
remaining = 1 - pv_coupons
s_T = 2 * ((1 + coupon) / remaining) ** (1/(2*T)) - 2
spot_rates[i] = s_T
# ---- Forward rates from spot rates ----------------------------- #
def forward_rate(s1, t1, s2, t2):
"""Continuously compounded forward rate between t1 and t2."""
return (s2 * t2 - s1 * t1) / (t2 - t1)
forward_rates = []
t_fwd = []
for i in range(1, len(maturities)):
f = forward_rate(spot_rates[i-1], maturities[i-1],
spot_rates[i], maturities[i])
forward_rates.append(f)
t_fwd.append((maturities[i-1] + maturities[i]) / 2)
# ---- Discount factors --------------------------------------- #
discount_factors = np.exp(-spot_rates * maturities)
df_result = pd.DataFrame({
'Maturity': maturities,
'Par Yield': par_yields * 100,
'Spot Rate': spot_rates * 100,
'Discount': discount_factors.round(6),
})
print(df_result.to_string(index=False))
# ---- Plot ---------------------------------------------------- #
fig, axes = plt.subplots(1, 2, figsize=(13, 5))
axes[0].plot(maturities, par_yields * 100, 'bo-', label='Par yields', linewidth=1.5, markersize=7)
axes[0].plot(maturities, spot_rates * 100, 'rs-', label='Spot rates', linewidth=1.5, markersize=7)
axes[0].plot(t_fwd, np.array(forward_rates) * 100, 'g^--', label='Forward rates', linewidth=1.5, markersize=6)
axes[0].set_xlabel("Maturity (years)"); axes[0].set_ylabel("Rate (%)")
axes[0].set_title("Yield Curve: Par, Spot, and Forward"); axes[0].legend(); axes[0].grid(alpha=0.3)
axes[1].plot(maturities, discount_factors, 'purple', linewidth=2.5)
axes[1].set_xlabel("Maturity (years)"); axes[1].set_ylabel("Discount Factor P(0,T)")
axes[1].set_title("Discount Factors"); axes[1].grid(alpha=0.3)
plt.tight_layout()
plt.savefig("yield_curve.png", dpi=150)
plt.show()
Step 3: Nelson-Siegel Yield Curve Fitting
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
# ------------------------------------------------------------------ #
# Fit Nelson-Siegel model to observed par yields
# ------------------------------------------------------------------ #
def nelson_siegel(m, beta0, beta1, beta2, lam):
"""
Nelson-Siegel yield curve.
m : maturity (years)
beta0 : long-run level
beta1 : slope
beta2 : curvature
lam : decay factor
"""
x = m / lam
factor1 = (1 - np.exp(-x)) / x
factor2 = factor1 - np.exp(-x)
return beta0 + beta1 * factor1 + beta2 * factor2
def svensson(m, b0, b1, b2, b3, l1, l2):
"""Svensson extension with second hump."""
x1 = m / l1; x2 = m / l2
f1 = (1 - np.exp(-x1)) / x1
f2_a = f1 - np.exp(-x1)
f3 = (1 - np.exp(-x2)) / x2 - np.exp(-x2)
return b0 + b1*f1 + b2*f2_a + b3*f3
# Observed yields (similar to US Treasury, 2023 inverted curve)
mats = np.array([0.25, 0.5, 1, 2, 3, 5, 7, 10, 20, 30])
yields_obs = np.array([5.30, 5.25, 5.10, 4.85, 4.65, 4.50, 4.40, 4.35, 4.50, 4.45]) / 100
# Fit Nelson-Siegel
def ns_loss(params):
b0, b1, b2, lam = params
if lam <= 0 or b0 <= 0:
return 1e10
y_fit = nelson_siegel(mats, b0, b1, b2, lam)
return np.sum((y_fit - yields_obs)**2)
result = minimize(ns_loss, x0=[0.045, -0.01, -0.02, 2.0],
method='Nelder-Mead', options={'maxiter': 10000, 'xatol': 1e-9})
b0, b1, b2, lam = result.x
print(f"Nelson-Siegel parameters:")
print(f" β₀ (level) = {b0*100:.4f}%")
print(f" β₁ (slope) = {b1*100:.4f}%")
print(f" β₂ (curvature)= {b2*100:.4f}%")
print(f" λ (decay) = {lam:.4f} years")
# In-sample fit quality
y_ns = nelson_siegel(mats, b0, b1, b2, lam)
rmse = np.sqrt(np.mean((y_ns - yields_obs)**2)) * 10000
print(f"RMSE: {rmse:.2f} bps")
# Plot
m_dense = np.linspace(0.25, 30, 200)
y_dense = nelson_siegel(m_dense, b0, b1, b2, lam)
fig, ax = plt.subplots(figsize=(10, 5))
ax.plot(mats, yields_obs * 100, 'ko', markersize=8, label='Observed', zorder=5)
ax.plot(m_dense, y_dense * 100, 'r-', linewidth=2.5, label=f'Nelson-Siegel (RMSE={rmse:.1f}bps)')
ax.set_xlabel("Maturity (years)"); ax.set_ylabel("Yield (%)")
ax.set_title("Nelson-Siegel Yield Curve Fitting")
ax.legend(); ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("nelson_siegel.png", dpi=150)
plt.show()
Advanced Usage
Credit Spread and CDS Pricing
import numpy as np
# ------------------------------------------------------------------ #
# CDS pricing: fair spread given hazard rate and recovery assumption
# ------------------------------------------------------------------ #
def cds_spread(lambda_h, R, T, dt=0.25):
"""
Compute par CDS spread (annualized).
Parameters
----------
lambda_h : float — constant hazard rate (default intensity)
R : float — recovery rate (e.g., 0.40)
T : float — CDS maturity (years)
dt : float — payment frequency (0.25 = quarterly)
"""
times = np.arange(dt, T + dt, dt)
# Survival probability: Q(t) = exp(-lambda_h * t)
Q = np.exp(-lambda_h * times)
# Risk-free discount: assume flat risk-free at 4%
r = 0.04
Z = np.exp(-r * times)
# Premium leg: sum of Q(t_i) * Z(t_i) * dt
pv_premium_leg = np.sum(Q * Z * dt)
# Protection leg: (1-R) * integral of Z(t) * lambda_h * Q(t)
# Approximation: sum over dt intervals
pv_protection_leg = (1 - R) * lambda_h * np.sum(Z * Q * dt)
return pv_protection_leg / pv_premium_leg
# Investment grade vs. high yield
for name, hazard, recovery in [
("AAA corporate", 0.0020, 0.40),
("BBB corporate", 0.0100, 0.40),
("BB high-yield", 0.0300, 0.40),
("B high-yield", 0.0600, 0.35),
("CCC distressed", 0.1500, 0.25),
]:
spread_5y = cds_spread(hazard, recovery, T=5) * 10000 # bps
print(f"{name:20s}: hazard={hazard*100:.2f}%/yr CDS spread={spread_5y:.0f} bps")
Troubleshooting
YTM solver fails (brentq ValueError)
Cause: Price is outside the valid range (too high/low for bounds).
Fix:
# Widen search bounds or check price
def safe_ytm(price, face, coupon, maturity):
try:
return brentq(lambda y: bond_price(face, coupon, y, maturity) - price,
-0.99, 50.0, xtol=1e-10)
except ValueError as e:
print(f"YTM search failed: {e}")
return np.nan
Nelson-Siegel fit doesn't converge
# Try multiple starting values
best = np.inf
best_params = None
for b0_0 in [0.03, 0.05, 0.07]:
for lam_0 in [1, 2, 5]:
res = minimize(ns_loss, x0=[b0_0, -0.01, 0.0, lam_0], method='Nelder-Mead')
if res.fun < best:
best = res.fun; best_params = res.x
Version Compatibility
| Package | Tested versions | Notes |
|---|---|---|
| scipy | 1.11, 1.12 | brentq and minimize stable |
| numpy | 1.24, 1.26 | No known issues |
External Resources
Official Documentation
Key Papers / Textbooks
- Fabozzi, F.J. (2012). Bond Markets, Analysis and Strategies (8th ed.). Pearson.
- Nelson, C.R. & Siegel, A.F. (1987). Parsimonious modeling of yield curves. Journal of Business.
Examples
Example 1: Immunization — Duration Matching
import numpy as np
from scipy.optimize import minimize
# Match duration of asset portfolio to liability duration
liability_pv = 1_000_000 # $1M liability in 7 years
liability_dur = 7.0 # years (zero-coupon liability)
# Available bonds: 2Y and 10Y
bond_2y = {'price': 960, 'duration': 1.90}
bond_10y = {'price': 900, 'duration': 7.50}
# Solve: w * dur_2 + (1-w) * dur_10 = liability_dur
# w * (dur_2 - dur_10) = liability_dur - dur_10
w = (liability_dur - bond_10y['duration']) / (bond_2y['duration'] - bond_10y['duration'])
print(f"Weight in 2Y bond: {w:.4f} ({w*100:.1f}%)")
print(f"Weight in 10Y bond: {1-w:.4f} ({(1-w)*100:.1f}%)")
pv_2y = w * liability_pv
pv_10y = (1-w) * liability_pv
n_2y = pv_2y / bond_2y['price']
n_10y = pv_10y / bond_10y['price']
print(f"Units of 2Y bond to buy: {n_2y:.0f}")
print(f"Units of 10Y bond to buy: {n_10y:.0f}")
Example 2: Key Rate Duration
import numpy as np
def key_rate_duration(bond, key_mats, bump=0.0001):
"""Compute key rate durations by bumping individual maturity points."""
krds = {}
for mat in key_mats:
# Bump yield at that maturity ±1bp
price_up = bond.price_at_yield(bond.ytm + bump) # simplified
price_down = bond.price_at_yield(bond.ytm - bump)
krd = -(price_up - price_down) / (2 * bump * bond.price)
krds[f"{mat}Y"] = krd
return krds
# For illustration with our Bond class
bond = Bond(1000, 0.05, 10, ytm=0.05)
# Single-factor KRD approximation (full KRD requires multi-factor spot curve)
print(f"Modified Duration (proxy for total KRD): {bond.modified_duration:.4f}")
print(f"Total DV01: ${bond.dv01:.4f} per $1000 face")
Last updated: 2026-03-17 | Maintainer: @xjtulyc Issues: GitHub Issues