"""bond-math-price-yield-duration-convexity.py — bond price, yield, duration, and convexity from first principles, applied to the Treasury constant-maturity curve as published on FRED. Input : datasets/bond-math-price-yield-duration-convexity.csv — the eleven nominal constant-maturity series DGS1MO DGS3MO DGS6MO DGS1 DGS2 DGS3 DGS5 DGS7 DGS10 DGS20 DGS30 (percent, daily, source Board of Governors H.15 via FRED), 2000-01-03 onward, retrieved 2026-09-11 through the public fredgraph CSV. Blank cells are FRED's "." (holidays; DGS1MO begins 2001-07-31). Re-download with --download. Method: A bond with face F, annual coupon rate c paid f times a year, and T years to maturity has n = f*T cash flows CF_k = F*c/f at times t_k = k/f, plus F at t_n. At a yield y compounded f times a year, P(y) = sum_k CF_k (1 + y/f)^(-f t_k) dP/dy = -sum_k CF_k t_k (1 + y/f)^(-f t_k - 1) d2P/dy2 = sum_k CF_k t_k (t_k + 1/f) (1 + y/f)^(-f t_k - 2) Macaulay duration D = sum_k t_k CF_k v_k / P (v_k = discount factor), modified duration Dm = -(dP/dy)/P = D/(1 + y/f), convexity C = (d2P/dy2)/P, DV01 = Dm * P * 0.0001. The second-order approximation of a price change for a yield shift dy is dP/P ~= -Dm dy + 0.5 C dy^2. Yield from price is found by bisection on P(y) - P_target. Every derivative is checked against a central finite difference before anything is printed; the script stops if the two disagree. The Treasury's constant-maturity curve is a par curve, so a bond whose coupon equals the quoted yield prices at 100 by construction; each maturity's yield is used as that maturity's coupon. Bill maturities (1, 3, 6 months) are treated as zero-coupon instruments discounted at the quoted investment-basis yield with semiannual compounding, so their Macaulay duration equals their maturity. Nothing here is a forecast, and nothing here is investment advice; the figures describe arithmetic on one day's published yields and on the published history of one series. Output: prints the article's tables; writes datasets/bond-math-price-yield-duration-convexity-shocks.csv (2-, 10-, 30-year par bonds repriced for yield shifts of -300 to +300 basis points, exact and approximate) and datasets/bond-math-price-yield-duration-convexity-history.csv (modified duration, convexity, and DV01 of a 10-year par bond at every DGS10 observation since 2000). With --figures, writes two PNG charts to content/images/bond-math-price-yield-duration-convexity/. Run : python code/bond-math-price-yield-duration-convexity.py [--download] [--figures] Needs : Python 3.13. Standard library only for every number. matplotlib 3.10.9 only for --figures. """ from __future__ import annotations import csv import io import math import os import sys import urllib.request ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) SLUG = "bond-math-price-yield-duration-convexity" DATA = os.path.join(ROOT, "datasets", SLUG + ".csv") SHOCKS = os.path.join(ROOT, "datasets", SLUG + "-shocks.csv") HISTORY = os.path.join(ROOT, "datasets", SLUG + "-history.csv") IMG_DIR = os.path.join(ROOT, "content", "images", SLUG) SERIES = ["DGS1MO", "DGS3MO", "DGS6MO", "DGS1", "DGS2", "DGS3", "DGS5", "DGS7", "DGS10", "DGS20", "DGS30"] YEARS = {"DGS1MO": 1 / 12, "DGS3MO": 0.25, "DGS6MO": 0.5, "DGS1": 1, "DGS2": 2, "DGS3": 3, "DGS5": 5, "DGS7": 7, "DGS10": 10, "DGS20": 20, "DGS30": 30} FACE = 100.0 FREQ = 2 # ------------------------------------------------------------------ the math --- def cash_flows(coupon: float, years: float, face: float = FACE, freq: int = FREQ) -> list[tuple[float, float]]: """(time in years, amount) for a level-coupon bond; a zero-coupon instrument is coupon=0.""" n = round(years * freq) if n == 0 or abs(n / freq - years) > 1e-9: # maturity is not a whole number of coupon periods: a single payment at T (a bill) return [(years, face)] flows = [(k / freq, face * coupon / freq) for k in range(1, n + 1)] t_n, cf_n = flows[-1] flows[-1] = (t_n, cf_n + face) return flows def price(coupon: float, ytm: float, years: float, face: float = FACE, freq: int = FREQ) -> float: return sum(cf * (1 + ytm / freq) ** (-freq * t) for t, cf in cash_flows(coupon, years, face, freq)) def dprice(coupon: float, ytm: float, years: float, face: float = FACE, freq: int = FREQ) -> float: return -sum(cf * t * (1 + ytm / freq) ** (-freq * t - 1) for t, cf in cash_flows(coupon, years, face, freq)) def d2price(coupon: float, ytm: float, years: float, face: float = FACE, freq: int = FREQ) -> float: return sum(cf * t * (t + 1 / freq) * (1 + ytm / freq) ** (-freq * t - 2) for t, cf in cash_flows(coupon, years, face, freq)) def macaulay(coupon: float, ytm: float, years: float, face: float = FACE, freq: int = FREQ) -> float: p = price(coupon, ytm, years, face, freq) return sum(t * cf * (1 + ytm / freq) ** (-freq * t) for t, cf in cash_flows(coupon, years, face, freq)) / p def modified(coupon: float, ytm: float, years: float, face: float = FACE, freq: int = FREQ) -> float: return -dprice(coupon, ytm, years, face, freq) / price(coupon, ytm, years, face, freq) def convexity(coupon: float, ytm: float, years: float, face: float = FACE, freq: int = FREQ) -> float: return d2price(coupon, ytm, years, face, freq) / price(coupon, ytm, years, face, freq) def yield_from_price(target: float, coupon: float, years: float, face: float = FACE, freq: int = FREQ, lo: float = -0.99, hi: float = 5.0, tol: float = 1e-12) -> float: """Bisection: P(y) is strictly decreasing in y, so the root is unique.""" f_lo = price(coupon, lo, years, face, freq) - target for _ in range(300): mid = 0.5 * (lo + hi) f_mid = price(coupon, mid, years, face, freq) - target if abs(f_mid) < tol: return mid if (f_lo > 0) == (f_mid > 0): lo, f_lo = mid, f_mid else: hi = mid return 0.5 * (lo + hi) def self_check() -> None: """The closed-form derivatives must agree with central finite differences, or nothing else is trusted.""" h = 1e-6 for coupon, ytm, years in ((0.0483, 0.0483, 10), (0.0528, 0.0528, 30), (0.02, 0.06, 7), (0.0, 0.04, 0.25)): p_up, p_dn = price(coupon, ytm + h, years), price(coupon, ytm - h, years) fd1 = (p_up - p_dn) / (2 * h) fd2 = (p_up - 2 * price(coupon, ytm, years) + p_dn) / (h * h) assert abs(fd1 - dprice(coupon, ytm, years)) < 1e-5, ("dP/dy", coupon, ytm, years) assert abs(fd2 - d2price(coupon, ytm, years)) / max(1.0, abs(fd2)) < 1e-3, ("d2P/dy2", coupon, ytm, years) assert abs(macaulay(coupon, ytm, years) / (1 + ytm / FREQ) - modified(coupon, ytm, years)) < 1e-12 assert abs(price(0.0483, 0.0483, 10) - 100.0) < 1e-9, "a par bond must price at 100" assert abs(yield_from_price(100.0, 0.0483, 10) - 0.0483) < 1e-10, "yield of a par bond is its coupon" # ------------------------------------------------------------------ the data --- def download() -> None: cols: dict[str, dict[str, str]] = {} for sid in SERIES: req = urllib.request.Request(f"https://fred.stlouisfed.org/graph/fredgraph.csv?id={sid}", headers={"User-Agent": "prismdatalab-research/1.0"}) with urllib.request.urlopen(req, timeout=60) as r: rows = list(csv.reader(io.StringIO(r.read().decode("utf-8")))) cols[sid] = {row[0]: ("" if row[1] == "." else row[1]) for row in rows[1:] if row} dates = sorted(d for d in set().union(*[set(c) for c in cols.values()]) if d >= "2000-01-01") with open(DATA, "w", newline="", encoding="utf-8") as fh: w = csv.writer(fh) w.writerow(["date"] + SERIES) for d in dates: w.writerow([d] + [cols[s].get(d, "") for s in SERIES]) print(f"downloaded {len(dates)} rows -> {DATA}") def load() -> list[dict[str, str]]: with open(DATA, encoding="utf-8") as fh: return list(csv.DictReader(fh)) def main() -> None: if "--download" in sys.argv: download() self_check() rows = load() latest = [r for r in rows if all(r[s] for s in SERIES)][-1] day = latest["date"] print(f"curve date {day} (dataset {rows[0]['date']} .. {rows[-1]['date']}, {len(rows)} rows)\n") # 1. the curve as a table of par bonds ----------------------------------------------------- print("Table 1: the constant-maturity curve as par instruments") print("series | years | yield % | price | Macaulay | modified | convexity | DV01 per 100") curve = {} for sid in SERIES: y = float(latest[sid]) / 100 T = YEARS[sid] coupon = 0.0 if T < 1 else y p = price(coupon, y, T) D, Dm, C = macaulay(coupon, y, T), modified(coupon, y, T), convexity(coupon, y, T) curve[sid] = (y, T, coupon, p, D, Dm, C) print(f"{sid} | {T:g} | {y*100:.2f} | {p:.4f} | {D:.3f} | {Dm:.3f} | {C:.2f} | {Dm*p*1e-4:.4f}") # 2. yield from price, the inverse problem ------------------------------------------------ y10 = curve["DGS10"][0] print("\nTable 2: yield solved from price, 10-year bond with coupon = the 2026 par yield") for target in (110.0, 105.0, 100.0, 95.0, 90.0): ysolved = yield_from_price(target, y10, 10) print(f"price {target:.2f} -> yield {ysolved*100:.4f} % (check: reprices to {price(y10, ysolved, 10):.6f})") # 3. shocks: exact vs first- and second-order --------------------------------------------- print("\nTable 3: exact repricing versus the duration and duration-plus-convexity approximations") print("bond | shift bp | exact price | duration only | dur+conv | err duration | err dur+conv") with open(SHOCKS, "w", newline="", encoding="utf-8") as fh: w = csv.writer(fh) w.writerow(["curve_date", "series", "years", "coupon_pct", "shift_bp", "yield_pct", "price_exact", "price_duration", "price_duration_convexity", "error_duration", "error_duration_convexity"]) for sid in ("DGS2", "DGS10", "DGS30"): y, T, coupon, p0, D, Dm, C = curve[sid] for bp in range(-300, 301, 25): dy = bp / 10000 exact = price(coupon, y + dy, T) p1 = p0 * (1 - Dm * dy) p2 = p0 * (1 - Dm * dy + 0.5 * C * dy * dy) w.writerow([day, sid, f"{T:g}", f"{coupon*100:.2f}", bp, f"{(y+dy)*100:.4f}", f"{exact:.6f}", f"{p1:.6f}", f"{p2:.6f}", f"{p1-exact:.6f}", f"{p2-exact:.6f}"]) if bp in (-300, -200, -100, -25, 25, 100, 200, 300): print(f"{sid} | {bp:+d} | {exact:.3f} | {p1:.3f} | {p2:.3f} | {p1-exact:+.3f} | {p2-exact:+.3f}") print(f"wrote {SHOCKS}") # 4. history: a 10-year par bond's risk numbers at every DGS10 print since 2000 ----------------- hist = [] with open(HISTORY, "w", newline="", encoding="utf-8") as fh: w = csv.writer(fh) w.writerow(["date", "DGS10_pct", "modified_duration_10y_par", "convexity_10y_par", "dv01_per_100"]) for r in rows: if not r["DGS10"]: continue y = float(r["DGS10"]) / 100 Dm, C = modified(y, y, 10), convexity(y, y, 10) hist.append((r["date"], y, Dm, C)) w.writerow([r["date"], f"{y*100:.2f}", f"{Dm:.4f}", f"{C:.3f}", f"{Dm*100*1e-4:.5f}"]) print(f"wrote {HISTORY} ({len(hist)} rows)") print("\nTable 4: modified duration of a 10-year par bond at selected DGS10 observations") lo = min(hist, key=lambda h: h[1]) hi = max(hist, key=lambda h: h[1]) for label, h in (("first observation", hist[0]), ("highest yield", hi), ("lowest yield", lo), ("latest observation", hist[-1])): print(f"{label} | {h[0]} | {h[1]*100:.2f} % | modified {h[2]:.3f} | convexity {h[3]:.2f} | DV01 {h[2]*1e-2:.4f}") yearly = {} for d, y, Dm, C in hist: yearly.setdefault(d[:4], []).append(Dm) print("\nTable 5: yearly mean modified duration of a 10-year par bond") for yr in sorted(yearly): print(f"{yr} | {sum(yearly[yr])/len(yearly[yr]):.2f}") if "--figures" in sys.argv: figures(curve, hist, day) # ---------------------------------------------------------------- the figures --- def figures(curve, hist, day) -> None: import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt os.makedirs(IMG_DIR, exist_ok=True) paper, ink, soft, grid = "#ffffff", "#1b2430", "#5e646b", "#e3e2de" c1, c2, c3 = "#2a78d6", "#eb6834", "#1baf7a" # fixed categorical order: 2-year, 10-year, 30-year plt.rcParams.update({"font.family": "DejaVu Sans", "font.size": 10, "axes.edgecolor": grid, "axes.labelcolor": ink, "xtick.color": soft, "ytick.color": soft, "text.color": ink, "axes.spines.top": False, "axes.spines.right": False}) # Figure 1: the 30-year par bond, exact price curve vs the two approximations, plus the errors fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11, 4.6), dpi=120, facecolor=paper) y, T, coupon, p0, D, Dm, C = curve["DGS30"] bps = list(range(-300, 301, 5)) exact = [price(coupon, y + b / 1e4, T) for b in bps] lin = [p0 * (1 - Dm * b / 1e4) for b in bps] quad = [p0 * (1 - Dm * b / 1e4 + 0.5 * C * (b / 1e4) ** 2) for b in bps] ax1.set_facecolor(paper) ax1.grid(axis="y", color=grid, linewidth=0.8) ax1.plot(bps, exact, color=c3, linewidth=2.2, label="exact repricing") ax1.plot(bps, quad, color=c2, linewidth=2, linestyle="--", label="duration + convexity") ax1.plot(bps, lin, color=c1, linewidth=2, linestyle=":", label="duration only") ax1.axhline(100, color=grid, linewidth=0.8) ax1.axvline(0, color=grid, linewidth=0.8) ax1.set_xlabel("parallel yield shift, basis points") ax1.set_ylabel("price per 100 face") ax1.set_title(f"30-year par bond, {y*100:.2f}% coupon ({day})", loc="left", fontsize=11) ax1.legend(frameon=False, loc="upper right") ax2.set_facecolor(paper) ax2.grid(axis="y", color=grid, linewidth=0.8) ax2.axhline(0, color=grid, linewidth=0.8) for sid, col, lab in (("DGS2", c1, "2-year"), ("DGS10", c2, "10-year"), ("DGS30", c3, "30-year")): y, T, coupon, p0, D, Dm, C = curve[sid] err1 = [p0 * (1 - Dm * b / 1e4) - price(coupon, y + b / 1e4, T) for b in bps] err2 = [p0 * (1 - Dm * b / 1e4 + 0.5 * C * (b / 1e4) ** 2) - price(coupon, y + b / 1e4, T) for b in bps] ax2.plot(bps, err1, color=col, linewidth=2, label=f"{lab}, duration only") ax2.plot(bps, err2, color=col, linewidth=2, linestyle="--", label=f"{lab}, with convexity") ax2.set_xlabel("parallel yield shift, basis points") ax2.set_ylabel("approximation minus exact, price points") ax2.set_title("Approximation error by maturity", loc="left", fontsize=11) ax2.legend(frameon=False, fontsize=8.5, ncol=2, loc="lower center") fig.tight_layout() out1 = os.path.join(IMG_DIR, "approximation-error.png") fig.savefig(out1, facecolor=paper) plt.close(fig) # Figure 2: yield and the resulting modified duration of a 10-year par bond, 2000-2026, stacked (one axis each) import datetime as dt dates = [dt.date.fromisoformat(h[0]) for h in hist] fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(11, 5.6), dpi=120, facecolor=paper, sharex=True) for ax in (ax1, ax2): ax.set_facecolor(paper) ax.grid(axis="y", color=grid, linewidth=0.8) ax1.plot(dates, [h[1] * 100 for h in hist], color=c1, linewidth=1.4) ax1.set_ylabel("DGS10, percent") ax1.set_title("10-year constant-maturity yield (FRED DGS10)", loc="left", fontsize=11) ax2.plot(dates, [h[2] for h in hist], color=c2, linewidth=1.4) ax2.set_ylabel("years") ax2.set_title("Modified duration of a 10-year par bond at that yield", loc="left", fontsize=11) fig.tight_layout() out2 = os.path.join(IMG_DIR, "duration-history.png") fig.savefig(out2, facecolor=paper) plt.close(fig) for out in (out1, out2): print(f"wrote {out} ({os.path.getsize(out) // 1024} KB)") if __name__ == "__main__": main()