"""nelson-siegel-yield-curve-fitting.py Nelson-Siegel and Svensson curves fitted to the eleven FRED constant-maturity Treasury yields, every trading day from 9 February 2006 (when the 30-year series came back) to the last day served. Prism Data Lab, 2026-10-13. Nothing here is investment advice. What it does 1. Reads the eleven daily constant-maturity series (DGS1MO ... DGS30, percent) and keeps the days on which all eleven have a value, from 2006-02-09. 2. Fits the Nelson-Siegel yield curve in the Diebold-Li form y(t) = b0 + b1 * (1 - exp(-L t)) / (L t) + b2 * ((1 - exp(-L t)) / (L t) - exp(-L t)), t in years, by ordinary least squares on each day, five ways: dl L fixed at Diebold and Li's 0.0609 per month (0.7308 per year) best L fixed at the one value that minimises the squared errors over every day (grid 0.050-3.000, 0.005) free L chosen afresh each day (grid 0.05-3.00, step 0.01), the nonlinear fit most papers describe svensson a second curvature term with its own decay L2 (both fixed over the sample, grid search), the extension the Federal Reserve's own curve (Gurkaynak, Sack and Wright, 2006) uses excl the dl curve fitted without the 1-month bill and the 20-year, the two maturities the misses point to, and judged at those two as out-of-sample points 3. Reports the errors (actual minus fitted, basis points) by model and maturity, the factors against Diebold and Li's data-based level, slope and curvature, the 20-year's miss before and after 20 May 2020, and one day, the last in the data, worked through. Inputs (one of) --download fetch the series from https://fred.stlouisfed.org/graph/fredgraph.csv?id=; robots.txt is read first and its Crawl-delay (at least one second) is kept between requests --raw-dir DIR read DIR/.csv, files saved from that endpoint earlier Outputs The tables, printed. With --csv, two files in --out (default: the datasets folder beside this script's folder): nelson-siegel-yield-curve-fitting.csv one row per model and maturity: mean error, MAE, RMSE, share above nelson-siegel-yield-curve-fitting-daily.csv one row per day: dl factors and RMSE, the free L and its RMSE, the svensson RMSE, and the 1-month and 20-year errors All eleven series are tagged "Public Domain: Citation Requested" on FRED (Board of Governors, H.15 release). With --figures, two charts in content/images/nelson-siegel-yield-curve-fitting/ (or --fig-out). Run python nelson-siegel-yield-curve-fitting.py --download --csv python nelson-siegel-yield-curve-fitting.py --raw-dir fred --csv --figures Standard library only for every number, Python 3.9 or later; matplotlib (3.10) only for --figures. About a minute. """ from __future__ import annotations import argparse import csv import io import math import operator import os import statistics import sys import time import urllib.request import urllib.robotparser SERIES = ["DGS1MO", "DGS3MO", "DGS6MO", "DGS1", "DGS2", "DGS3", "DGS5", "DGS7", "DGS10", "DGS20", "DGS30"] TAU = [1 / 12, 0.25, 0.5, 1, 2, 3, 5, 7, 10, 20, 30] # years LABEL = ["1M", "3M", "6M", "1Y", "2Y", "3Y", "5Y", "7Y", "10Y", "20Y", "30Y"] START = "2006-02-09" # FRED: the 30-year series was discontinued 2002-02-18 and reintroduced 2006-02-09; # H.15 footnote 9: the values between are Treasury estimates from the 20-year LAM_DL = 0.0609 * 12 # Diebold and Li (2006): 0.0609 with maturity in months TWENTY_BACK = "2020-05-20" # Treasury: the 20-year bond was reintroduced on May 20, 2020 GRID_BEST = [round(0.05 + 0.005 * i, 3) for i in range(591)] # 0.050 .. 3.000 GRID_FREE = [round(0.05 + 0.01 * i, 2) for i in range(296)] # 0.05 .. 3.00 GRID_L1 = [round(0.05 * i, 2) for i in range(1, 61)] # 0.05 .. 3.00 GRID_L2 = [round(0.01 * i, 2) for i in range(1, 101)] # 0.01 .. 1.00 USER_AGENT = "nelson-siegel-fit/1.0 (put your own contact here)" SLUG = "nelson-siegel-yield-curve-fitting" EXCL = [0, 9] # indexes of DGS1MO and DGS20 # ----------------------------------------------------------------- data --- def parse_fredgraph(text: str) -> dict[str, float]: rows = list(csv.reader(io.StringIO(text))) out = {} for r in rows[1:]: if r and r[1].strip() not in ("", "."): out[r[0]] = float(r[1]) return out def load_all(raw_dir: str | None) -> dict[str, dict[str, float]]: data = {} if raw_dir: for sid in SERIES: with open(os.path.join(raw_dir, sid + ".csv"), encoding="utf-8") as f: data[sid] = parse_fredgraph(f.read()) return data headers = {"User-Agent": USER_AGENT, "Accept": "*/*"} # FRED stalls on a request with no Accept header req = urllib.request.Request("https://fred.stlouisfed.org/robots.txt", headers=headers) with urllib.request.urlopen(req, timeout=60) as r: rules = r.read().decode("utf-8", "replace") rp = urllib.robotparser.RobotFileParser() rp.parse(rules.splitlines()) delay = max(1.0, float(rp.crawl_delay(USER_AGENT) or 0)) for sid in SERIES: url = "https://fred.stlouisfed.org/graph/fredgraph.csv?id=" + sid if not rp.can_fetch(USER_AGENT, url): sys.exit(f"robots.txt disallows {url}; stopping") with urllib.request.urlopen(urllib.request.Request(url, headers=headers), timeout=90) as r: data[sid] = parse_fredgraph(r.read().decode("utf-8")) time.sleep(delay) return data def panel(data) -> tuple[list[str], list[list[float]]]: days = sorted(d for d in set().union(*[set(v) for v in data.values()]) if d >= START and all(d in data[s] for s in SERIES)) return days, [[data[s][d] for s in SERIES] for d in days] # ---------------------------------------------------------- linear algebra --- def ns_row(t: float, lam: float) -> list[float]: x = lam * t s = (1 - math.exp(-x)) / x return [1.0, s, s - math.exp(-x)] def sv_row(t: float, l1: float, l2: float) -> list[float]: x2 = l2 * t s2 = (1 - math.exp(-x2)) / x2 return ns_row(t, l1) + [s2 - math.exp(-x2)] def solve(a: list[list[float]], b: list[list[float]]) -> list[list[float]]: """Solve A X = B (A square) by Gaussian elimination with partial pivoting.""" n = len(a) m = [row[:] + brow[:] for row, brow in zip(a, b)] for c in range(n): p = max(range(c, n), key=lambda r: abs(m[r][c])) m[c], m[p] = m[p], m[c] for r in range(n): if r != c: f = m[r][c] / m[c][c] m[r] = [x - f * y for x, y in zip(m[r], m[c])] return [[x / m[r][r] for x in m[r][n:]] for r in range(n)] def hat(x: list[list[float]]) -> tuple[list[list[float]], list[list[float]]]: """For a design X (n x k): the OLS operator (X'X)^-1 X' (k x n) and the projection X (X'X)^-1 X' (n x n).""" k = len(x[0]) xtx = [[sum(r[i] * r[j] for r in x) for j in range(k)] for i in range(k)] xt = [[r[i] for r in x] for i in range(k)] op = solve(xtx, xt) proj = [[sum(x[i][c] * op[c][j] for c in range(k)) for j in range(len(x))] for i in range(len(x))] return op, proj def sse_total(proj, cross) -> float: """Sum over days of y'(I - P)y = trace(M) - sum_ij P_ij M_ij, with M = sum of y y' over days.""" n = len(proj) return sum(cross[i][i] for i in range(n)) - sum(proj[i][j] * cross[i][j] for i in range(n) for j in range(n)) def fit_days(Y, x): op, _ = hat(x) betas, resid = [], [] for y in Y: b = [sum(o * v for o, v in zip(row, y)) for row in op] fitted = [sum(bb * xx for bb, xx in zip(b, xr)) for xr in x] betas.append(b) resid.append([v - f for v, f in zip(y, fitted)]) return betas, resid def curvature_peak() -> float: """x maximising (1 - e^-x)/x - e^-x, by bisection on its derivative.""" d = lambda x: (math.exp(-x) * (x * x + x + 1) - 1) / (x * x) # noqa: E731 derivative of the loading lo, hi = 0.5, 5.0 for _ in range(200): mid = (lo + hi) / 2 lo, hi = (mid, hi) if d(mid) > 0 else (lo, mid) return (lo + hi) / 2 # ----------------------------------------------------------- statistics --- def rmse(e): return math.sqrt(sum(v * v for v in e) / len(e)) def corr(a, b): ma, mb = statistics.fmean(a), statistics.fmean(b) sab = sum((x - ma) * (y - mb) for x, y in zip(a, b)) return sab / math.sqrt(sum((x - ma) ** 2 for x in a) * sum((y - mb) ** 2 for y in b)) def pct(xs, q): s = sorted(xs) pos = (len(s) - 1) * q lo = int(pos) return s[lo] + (s[min(lo + 1, len(s) - 1)] - s[lo]) * (pos - lo) def by_maturity(resid, cols=range(11)): out = {} for j in cols: e = [100 * r[j] for r in resid] out[j] = {"mean": statistics.fmean(e), "mae": statistics.fmean(abs(v) for v in e), "rmse": rmse(e), "above": sum(1 for v in e if v > 0) / len(e)} return out # ----------------------------------------------------------------- run --- def run(data): days, Y = panel(data) n = len(days) cross = [[sum(y[i] * y[j] for y in Y) for j in range(11)] for i in range(11)] res = {"days": days, "Y": Y} # dl: Diebold-Li's fixed decay x_dl = [ns_row(t, LAM_DL) for t in TAU] res["dl_beta"], res["dl_resid"] = fit_days(Y, x_dl) # best: one decay for the whole sample, by total squared error sse = [(sse_total(hat([ns_row(t, lam) for t in TAU])[1], cross), lam) for lam in GRID_BEST] res["best_lam"] = min(sse)[1] res["best_beta"], res["best_resid"] = fit_days(Y, [ns_row(t, res["best_lam"]) for t in TAU]) # free: a decay per day; y'Py from the 66 distinct products y_i y_j pairs = [(i, j) for i in range(11) for j in range(i, 11)] pw = [] for lam in GRID_FREE: proj = hat([ns_row(t, lam) for t in TAU])[1] pw.append([proj[i][j] * (1 if i == j else 2) for i, j in pairs]) free_lam, free_rmse = [], [] for y in Y: yy = [y[i] * y[j] for i, j in pairs] tot = sum(v * v for v in y) best = max(range(len(GRID_FREE)), key=lambda k: sum(map(operator.mul, pw[k], yy))) free_lam.append(GRID_FREE[best]) free_rmse.append(100 * math.sqrt(max(tot - sum(map(operator.mul, pw[best], yy)), 0.0) / 11)) res["free_lam"], res["free_rmse"] = free_lam, free_rmse # svensson: two decays, both fixed over the sample (L2 below L1, so the second hump sits further out) sv = [] for l1 in GRID_L1: for l2 in GRID_L2: if round(l2 * 100) <= round(l1 * 100) - 5: # in hundredths, so no float edge drops a pair sv.append((sse_total(hat([sv_row(t, l1, l2) for t in TAU])[1], cross), l1, l2)) _, res["sv_l1"], res["sv_l2"] = min(sv) res["sv_candidates"] = len(sv) res["sv_beta"], res["sv_resid"] = fit_days(Y, [sv_row(t, res["sv_l1"], res["sv_l2"]) for t in TAU]) # excl: dl fitted on nine maturities, judged at the 1-month and the 20-year keep = [j for j in range(11) if j not in EXCL] x9 = [ns_row(TAU[j], LAM_DL) for j in keep] b9, _ = fit_days([[y[j] for j in keep] for y in Y], x9) res["excl_resid"] = [[y[j] - sum(bb * xx for bb, xx in zip(b, ns_row(TAU[j], LAM_DL))) for j in range(11)] for y, b in zip(Y, b9)] res["excl_fit_rmse"] = [100 * rmse([r[j] for j in keep]) for r in res["excl_resid"]] res["n"] = n return res def report(res): days, Y, n = res["days"], res["Y"], res["n"] print(f"Sample: {n:,} trading days with all eleven yields, {days[0]} to {days[-1]}") xpk = curvature_peak() print(f"Curvature loading peaks at L*t = {xpk:.4f}: Diebold-Li's L puts it at {xpk / 0.0609:.2f} months; " f"the best fixed L ({res['best_lam']:.3f} per year) at {xpk / res['best_lam']:.2f} years") models = [("dl", f"Nelson-Siegel, L = {LAM_DL:.4f}/yr (Diebold-Li)", res["dl_resid"]), ("best", f"Nelson-Siegel, L = {res['best_lam']:.3f}/yr (best fixed)", res["best_resid"]), ("svensson", f"Svensson, L1 = {res['sv_l1']:.2f}, L2 = {res['sv_l2']:.2f}/yr", res["sv_resid"])] daily = {} for key, name, resid in models: r = [100 * rmse(e) for e in resid] daily[key] = r print(f"\n{name}: daily RMSE median {statistics.median(r):.2f} bp, mean {statistics.fmean(r):.2f}, " f"90th pct {pct(r, 0.9):.2f}, max {max(r):.2f} on {days[r.index(max(r))]}") bm = by_maturity(resid) print(" maturity mean MAE RMSE above") for j in range(11): s = bm[j] print(f" {LABEL[j]:>4} {s['mean']:+7.2f} {s['mae']:6.2f} {s['rmse']:6.2f} {100 * s['above']:5.1f}%") print(f"\nSvensson grid: {res['sv_candidates']:,} decay pairs searched; humps at {xpk / res['sv_l1']:.2f} and " f"{xpk / res['sv_l2']:.2f} years") fr = res["free_rmse"] edge = sum(1 for v in res["free_lam"] if v in (GRID_FREE[0], GRID_FREE[-1])) steps = [abs(a - b) for a, b in zip(res["free_lam"][1:], res["free_lam"])] print(f"\nFree L per day: RMSE median {statistics.median(fr):.2f} bp, 90th pct {pct(fr, 0.9):.2f}; L 10th/50th/90th pct " f"{pct(res['free_lam'], .1):.2f} / {pct(res['free_lam'], .5):.2f} / {pct(res['free_lam'], .9):.2f}; at a grid " f"edge on {edge} days; median day-to-day move {statistics.median(steps):.2f}, moves of 0.5 or more on " f"{sum(1 for s in steps if s >= 0.5)} days") ex = res["excl_resid"] print(f"\nNine-maturity fit (no 1M, no 20Y): in-sample RMSE median {statistics.median(res['excl_fit_rmse']):.2f} bp") for j in EXCL: e = [100 * r[j] for r in ex] print(f" {LABEL[j]} out-of-sample: mean {statistics.fmean(e):+.2f} bp, above the curve on " f"{100 * sum(1 for v in e if v > 0) / n:.1f}% of days") for lab, sel in (("before", lambda d: d < TWENTY_BACK), ("from", lambda d: d >= TWENTY_BACK)): idx = [i for i, d in enumerate(days) if sel(d)] e = [100 * ex[i][9] for i in idx] inv = sum(1 for i in idx if Y[i][9] > Y[i][10]) print(f" 20Y {lab} {TWENTY_BACK}: {len(idx):,} days, mean out-of-sample miss {statistics.fmean(e):+.2f} bp, " f"above the curve on {100 * sum(1 for v in e if v > 0) / len(e):.1f}% of days, smallest {min(e):+.2f}; " f"20Y above 30Y on {inv:,} days ({100 * inv / len(idx):.1f}%)") e20 = [100 * r[9] for r in ex] k = e20.index(max(e20)) print(f" 20Y largest out-of-sample miss {e20[k]:+.2f} bp on {days[k]}; by calendar year (mean, bp):") years = sorted({d[:4] for d in days}) print(" " + " ".join(f"{y} {statistics.fmean(e20[i] for i, d in enumerate(days) if d[:4] == y):+.1f}" for y in years)) b = res["dl_beta"] lvl = [y[8] for y in Y] slope = [y[8] - y[1] for y in Y] curv = [2 * y[4] - y[1] - y[8] for y in Y] print(f"\nFactors vs Diebold-Li's data-based measures: corr(b0, 10Y) {corr([x[0] for x in b], lvl):.3f}; " f"corr(b1, 10Y-3M) {corr([x[1] for x in b], slope):.3f}; corr(b2, 2*2Y-3M-10Y) " f"{corr([x[2] for x in b], curv):.3f}; corr(b0 + b1, 3M) {corr([x[0] + x[1] for x in b], [y[1] for y in Y]):.3f}") i = n - 1 xs = [ns_row(t, LAM_DL) for t in TAU] print(f"\nWorked example, {days[i]} (Diebold-Li L): b0 {b[i][0]:.4f}, b1 {b[i][1]:.4f}, b2 {b[i][2]:.4f}; " f"RMSE {daily['dl'][i]:.2f} bp") print(" maturity actual slope-load curv-load fitted error(bp)") for j in range(11): f = sum(bb * xx for bb, xx in zip(b[i], xs[j])) print(f" {LABEL[j]:>4} {Y[i][j]:7.2f} {xs[j][1]:10.4f} {xs[j][2]:10.4f} {f:7.3f} {100 * (Y[i][j] - f):+8.1f}") return daily def write_csv(res, daily, out): os.makedirs(out, exist_ok=True) p1 = os.path.join(out, SLUG + ".csv") with open(p1, "w", newline="", encoding="utf-8") as f: w = csv.writer(f) w.writerow(["model", "decay_per_year", "maturity", "series", "days", "mean_error_bp", "mae_bp", "rmse_bp", "share_above"]) for key, decay, resid, cols in (("dl", f"{LAM_DL:.4f}", res["dl_resid"], range(11)), ("best", f"{res['best_lam']:.3f}", res["best_resid"], range(11)), ("svensson", f"{res['sv_l1']:.2f};{res['sv_l2']:.2f}", res["sv_resid"], range(11)), ("excl", f"{LAM_DL:.4f}", res["excl_resid"], range(11))): for j, s in by_maturity(resid, cols).items(): w.writerow([key, decay, LABEL[j], SERIES[j], res["n"], f"{s['mean']:.2f}", f"{s['mae']:.2f}", f"{s['rmse']:.2f}", f"{s['above']:.4f}"]) p2 = os.path.join(out, SLUG + "-daily.csv") with open(p2, "w", newline="", encoding="utf-8") as f: w = csv.writer(f) w.writerow(["date", "b0", "b1", "b2", "rmse_dl_bp", "free_decay", "rmse_free_bp", "rmse_svensson_bp", "error_1m_bp", "error_20y_bp", "error_20y_excl_bp"]) for i, d in enumerate(res["days"]): b = res["dl_beta"][i] w.writerow([d, f"{b[0]:.4f}", f"{b[1]:.4f}", f"{b[2]:.4f}", f"{daily['dl'][i]:.2f}", f"{res['free_lam'][i]:.2f}", f"{res['free_rmse'][i]:.2f}", f"{daily['svensson'][i]:.2f}", f"{100 * res['dl_resid'][i][0]:.2f}", f"{100 * res['dl_resid'][i][9]:.2f}", f"{100 * res['excl_resid'][i][9]:.2f}"]) print(f"wrote {p1}\nwrote {p2}") def figures(res, out_dir): import matplotlib matplotlib.use("Agg") import matplotlib.dates as mdates import matplotlib.pyplot as plt from datetime import date paper, ink, muted, grid, axis = "#ffffff", "#101820", "#56616b", "#e3e6e8", "#9aa3ab" ramp = ["#8fbcc2", "#4a8a94", "#0b5563"] # the site's validated one-hue ramp, light to dark minus = lambda s: s.replace("-", "−") # noqa: E731 plt.rcParams.update({"font.family": "DejaVu Sans", "font.size": 10, "text.color": ink, "xtick.color": muted, "ytick.color": muted}) os.makedirs(out_dir, exist_ok=True) # 1: mean error by maturity, three curves fig, ax = plt.subplots(figsize=(8.8, 4.6), dpi=150, facecolor=paper) fig.subplots_adjust(left=0.085, right=0.985, top=0.80, bottom=0.11) sets = [("Nelson-Siegel, Diebold-Li decay", res["dl_resid"]), (f"Nelson-Siegel, best fixed decay ({res['best_lam']:.2f}/yr)", res["best_resid"]), ("Svensson (second hump)", res["sv_resid"])] wbar = 0.26 for k, (name, resid) in enumerate(sets): bm = by_maturity(resid) ax.bar([j + (k - 1) * wbar for j in range(11)], [bm[j]["mean"] for j in range(11)], width=wbar - 0.03, color=ramp[k], label=name, linewidth=0, zorder=3) ax.axhline(0, color=axis, linewidth=1.0, zorder=2) ax.yaxis.grid(True, color=grid, linewidth=0.8, zorder=0) for s in ax.spines.values(): s.set_visible(False) ax.set_xticks(range(11), LABEL) ax.tick_params(length=0, labelsize=9.5) ax.yaxis.set_major_formatter(matplotlib.ticker.FuncFormatter(lambda v, _: minus(f"{v:+g}") if v else "0")) ax.set_ylabel("basis points", color=muted) fig.text(0.085, 0.95, f"Average miss by maturity, actual minus fitted, {res['n']:,} trading days", fontsize=12.5, ha="left", va="center") fig.legend(loc="center left", bbox_to_anchor=(0.078, 0.875), ncol=3, frameon=False, fontsize=9, handlelength=1.1, columnspacing=1.4) p = os.path.join(out_dir, "average-miss-by-maturity.png") fig.savefig(p, facecolor=paper) plt.close(fig) print(f"wrote {p} ({os.path.getsize(p) / 1024:.0f} KB)") # 2: the 20-year's out-of-sample miss over time (nine-maturity fit), 21-day mean days = [date.fromisoformat(d) for d in res["days"]] e = [100 * r[9] for r in res["excl_resid"]] sm = [statistics.fmean(e[max(0, i - 20):i + 1]) for i in range(len(e))] fig, ax = plt.subplots(figsize=(8.8, 4.2), dpi=150, facecolor=paper) fig.subplots_adjust(left=0.085, right=0.985, top=0.86, bottom=0.10) ax.plot(days, e, color=ramp[0], linewidth=0.6, zorder=2) ax.plot(days, sm, color=ramp[2], linewidth=2.0, zorder=3) ax.axhline(0, color=axis, linewidth=1.0, zorder=1) back = date.fromisoformat(TWENTY_BACK) ax.axvline(back, color=muted, linewidth=1.0, linestyle=(0, (3, 3)), zorder=1) top = max(e) ax.annotate("20-year bond reissued, 20 May 2020", (back, top), xytext=(-6, 0), textcoords="offset points", ha="right", va="top", fontsize=9, color=muted) ax.yaxis.grid(True, color=grid, linewidth=0.8, zorder=0) for s in ax.spines.values(): s.set_visible(False) ax.tick_params(length=0, labelsize=9.5) ax.xaxis.set_major_locator(mdates.YearLocator(2)) ax.xaxis.set_major_formatter(mdates.DateFormatter("%Y")) ax.yaxis.set_major_formatter(matplotlib.ticker.FuncFormatter(lambda v, _: minus(f"{v:+g}") if v else "0")) ax.set_ylabel("basis points", color=muted) fig.text(0.085, 0.94, "The 20-year against a curve fitted without it: daily (pale) and 21-day mean (dark)", fontsize=12.5, ha="left", va="center") p = os.path.join(out_dir, "twenty-year-out-of-sample-miss.png") fig.savefig(p, facecolor=paper) plt.close(fig) print(f"wrote {p} ({os.path.getsize(p) / 1024:.0f} KB)") def main(): ap = argparse.ArgumentParser(description=__doc__.split("\n\n")[0]) src = ap.add_mutually_exclusive_group(required=True) src.add_argument("--download", action="store_true", help="fetch the eleven series from FRED") src.add_argument("--raw-dir", help="read /.csv saved from FRED's fredgraph.csv endpoint") ap.add_argument("--csv", action="store_true", help="write the two CSV files") here = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) ap.add_argument("--out", default=os.path.join(here, "datasets")) ap.add_argument("--figures", action="store_true", help="draw the two charts (needs matplotlib)") ap.add_argument("--fig-out", default=os.path.join(here, "content", "images", SLUG)) a = ap.parse_args() sys.stdout.reconfigure(encoding="utf-8") res = run(load_all(a.raw_dir)) daily = report(res) if a.csv: write_csv(res, daily, a.out) if a.figures: figures(res, a.fig_out) if __name__ == "__main__": main()