The Nelson-Siegel curve¶

The yield at a tenor from the curve's three factors, y(t) = b0 + b1 L + b2 (L - e^-x), with x = t / tau and the slope's loading L = (1 - e^-x) / x: the long end tends to the level b0, the short end to the level plus the slope, and the curvature bends the middle.

Below x = 1/2, where 1 - e^-x in binary32 would lose the loading's digits, L is its series, the sum of (-x)^k / (k + 1)!, to the term whose next is below the float's rounding there. Source: Nelson & Siegel, Parsimonious Modeling of Yield Curves, J. Business 60 (1987); Diebold & Li, J. Econometrics 130 (2006). Checked against a reference in numpy, which has its own tests.

Source

Nelson & Siegel, Parsimonious Modeling of Yield Curves, J. Business 60 (1987); Diebold & Li, J. Econometrics 130 (2006).

Proof

curve, limit.

nelson_siegel.trait

namespace lusine.finance;

float<year^-1> nelson_siegel(float<year> tenor, float<year^-1> level, float<year^-1> slope, float<year^-1> curvature,
                             float<year> tau)
{
  let x = tenor / tau;
  let e = exp(-x);
  let series = 1.0 - x / 2.0 * (1.0 - x / 3.0 * (1.0 - x / 4.0 * (1.0 - x / 5.0 * (1.0 - x / 6.0 * (1.0 - x / 7.0
               * (1.0 - x / 8.0 * (1.0 - x / 9.0)))))));
  let load = x < 0.5 ? series : (1.0 - e) / x;
  return level + slope * load + curvature * (load - e);
}

The yield at a tenor from the curve's three factors, level, slope and curvature, y(t) = b0 + b1 (1 - e^-x) / x + b2 ((1 - e^-x) / x - e^-x), x = t / tau.

In [1]:

import os, sys
root = os.getcwd()
while not os.path.isfile(os.path.join(root, 'docsite.json')):
    root = os.path.dirname(root)
sys.path.insert(0, os.path.join(root, 'tools'))
import numpy as np
import matplotlib.pyplot as plt
from lusine import trait_run
EPS = float(np.finfo(np.float32).eps)
YEAR = 31557600.0      # lusine's year row, the Julian year in seconds: the host stores and reads in the base units

import os, sys root = os.getcwd() while not os.path.isfile(os.path.join(root, 'docsite.json')): root = os.path.dirname(root) sys.path.insert(0, os.path.join(root, 'tools')) import numpy as np import matplotlib.pyplot as plt from lusine import trait_run EPS = float(np.finfo(np.float32).eps) YEAR = 31557600.0 # lusine's year row, the Julian year in seconds: the host stores and reads in the base units

The kernel¶

The law as the library states it, and as this page runs it.

In [2]:

print(trait_run.source('finance/nelson_siegel.trait'))

print(trait_run.source('finance/nelson_siegel.trait'))

// lusine.finance.nelson_siegel: the yield at a tenor from the curve's three factors,
// y(t) = b0 + b1 L + b2 (L - e^-x), with x = t / tau and the slope's loading L = (1 - e^-x) / x: the long end tends to
// the level b0, the short end to the level plus the slope, and the curvature bends the middle. Below x = 1/2, where
// 1 - e^-x in binary32 would lose the loading's digits, L is its series, the sum of (-x)^k / (k + 1)!, to the term
// whose next is below the float's rounding there.
// Source: Nelson & Siegel, Parsimonious Modeling of Yield Curves, J. Business 60 (1987); Diebold & Li,
// J. Econometrics 130 (2006).
// Checked against a reference in numpy, which has its own tests.

namespace lusine.finance;

float<year^-1> nelson_siegel(float<year> tenor, float<year^-1> level, float<year^-1> slope, float<year^-1> curvature,
                             float<year> tau)
{
  let x = tenor / tau;
  let e = exp(-x);
  let series = 1.0 - x / 2.0 * (1.0 - x / 3.0 * (1.0 - x / 4.0 * (1.0 - x / 5.0 * (1.0 - x / 6.0 * (1.0 - x / 7.0
               * (1.0 - x / 8.0 * (1.0 - x / 9.0)))))));
  let load = x < 0.5 ? series : (1.0 - e) / x;
  return level + slope * load + curvature * (load - e);
}

The program¶

Each section runs the kernel over a lattice of items, one a case: the host sets every parameter, and main() binds the kernel over the items. The host speaks the base units, so a time in years is stored as t * YEAR seconds and an annual rate as r / YEAR, and a result in years is read back the same way; inside, the law runs on those, in binary32. curve below is that call; closed is the independent account, the same law in double precision with the loading from numpy's expm1. Each yield is within 4 ulps of its largest term.

In [3]:

print(trait_run.law_program('finance/nelson_siegel.trait', 'nelson_siegel'))

def curve(tenor, level, slope, curvature, tau=1.5):
    """The kernel's yield at each tenor, the tenors and tau in years, the factors and the yield a year."""
    return trait_run.law('finance/nelson_siegel.trait', 'nelson_siegel', tenor=np.asarray(tenor) * YEAR,
                         level=np.asarray(level) / YEAR, slope=np.asarray(slope) / YEAR,
                         curvature=np.asarray(curvature) / YEAR, tau=np.asarray(tau) * YEAR) * YEAR


def closed(tenor, level, slope, curvature, tau=1.5):
    """The independent account, in double precision, the loading by numpy's expm1."""
    x = np.asarray(tenor, dtype=float) / tau
    with np.errstate(invalid='ignore', divide='ignore'):
        load = np.where(x > 0, -np.expm1(-x) / x, 1.0)
    return level + slope * load + curvature * (load - np.exp(-x))

print(trait_run.law_program('finance/nelson_siegel.trait', 'nelson_siegel')) def curve(tenor, level, slope, curvature, tau=1.5): """The kernel's yield at each tenor, the tenors and tau in years, the factors and the yield a year.""" return trait_run.law('finance/nelson_siegel.trait', 'nelson_siegel', tenor=np.asarray(tenor) * YEAR, level=np.asarray(level) / YEAR, slope=np.asarray(slope) / YEAR, curvature=np.asarray(curvature) / YEAR, tau=np.asarray(tau) * YEAR) * YEAR def closed(tenor, level, slope, curvature, tau=1.5): """The independent account, in double precision, the loading by numpy's expm1.""" x = np.asarray(tenor, dtype=float) / tau with np.errstate(invalid='ignore', divide='ignore'): load = np.where(x > 0, -np.expm1(-x) / x, 1.0) return level + slope * load + curvature * (load - np.exp(-x))

namespace probe;

lusine.version version = { { number = 1 } };

struct item {
  host float<year> tenor;
  host float<year^-1> level;
  host float<year^-1> slope;
  host float<year^-1> curvature;
  host float<year> tau;
  float<year^-1> result_0;
};
lattice<item, 1> items(spacing = 1 m, boundary = wall);

void main() { items.result_0 = lusine.finance.nelson_siegel(items.tenor, items.level, items.slope, items.curvature, items.tau); }

Curve: normal, inverted and humped¶

Three factor sets, 300 tenors each in one lattice: a rising curve, a falling one, and a hump whose peak is inside the tenors.

In [4]:

t = np.linspace(0.05, 30, 300)
shapes = {'normal': (0.045, -0.02, 0.0), 'inverted': (0.03, 0.02, 0.0), 'humped': (0.04, -0.01, 0.04)}
f = np.repeat(np.array(list(shapes.values())), len(t), axis=0)
y = curve(np.tile(t, 3), f[:, 0], f[:, 1], f[:, 2]).reshape(3, len(t))
fig, ax = plt.subplots(figsize=(7, 3.5))
for (name, (b0, b1, b2)), yy in zip(shapes.items(), y):
    assert np.all(np.abs(yy - closed(t, b0, b1, b2)) <= 4 * EPS * (abs(b0) + abs(b1) + abs(b2)))
    ax.plot(t, yy * 100, label=name)
assert np.all(np.diff(y[0]) > 0) and np.all(np.diff(y[1]) < 0) and 0 < np.argmax(y[2]) < len(t) - 1
ax.set_xlabel('tenor (years)')
ax.set_ylabel('yield (%)')
ax.legend()
plt.show()

t = np.linspace(0.05, 30, 300) shapes = {'normal': (0.045, -0.02, 0.0), 'inverted': (0.03, 0.02, 0.0), 'humped': (0.04, -0.01, 0.04)} f = np.repeat(np.array(list(shapes.values())), len(t), axis=0) y = curve(np.tile(t, 3), f[:, 0], f[:, 1], f[:, 2]).reshape(3, len(t)) fig, ax = plt.subplots(figsize=(7, 3.5)) for (name, (b0, b1, b2)), yy in zip(shapes.items(), y): assert np.all(np.abs(yy - closed(t, b0, b1, b2)) <= 4 * EPS * (abs(b0) + abs(b1) + abs(b2))) ax.plot(t, yy * 100, label=name) assert np.all(np.diff(y[0]) > 0) and np.all(np.diff(y[1]) < 0) and 0 < np.argmax(y[2]) < len(t) - 1 ax.set_xlabel('tenor (years)') ax.set_ylabel('yield (%)') ax.legend() plt.show()

Limit: the short end is level + slope, the long end the level¶

As the tenor falls to nothing the curve tends to level + slope, and at no tenor it is that; as the tenor grows, the curve tends to the level, the slope's and the curvature's loadings both tau / t.

In [5]:

lv, sl, cu = 0.04, -0.015, 0.02
y = curve([0.0, 1e-8, 1e5], lv, sl, cu)
assert np.all(np.abs(y[:2] - (lv + sl)) <= 2 * EPS * lv), y[:2]
assert abs(y[2] - (lv + (sl + cu) * 1.5 / 1e5)) <= 2 * EPS * lv, y[2]
print('short end', y[:2], ' long end', y[2])

lv, sl, cu = 0.04, -0.015, 0.02 y = curve([0.0, 1e-8, 1e5], lv, sl, cu) assert np.all(np.abs(y[:2] - (lv + sl)) <= 2 * EPS * lv), y[:2] assert abs(y[2] - (lv + (sl + cu) * 1.5 / 1e5)) <= 2 * EPS * lv, y[2] print('short end', y[:2], ' long end', y[2])

short end [0.025 0.025]  long end 0.04000007764775049

Invariant: one curve across the series' seam, and linear in its factors¶

Below x = 1/2 the loading is its series, above it (1 - e^-x) / x: across the seam the curve is continuous to binary32's rounding, each side within the account's tolerance. Doubling every factor doubles the yield, bit for bit.

In [6]:

lv, sl, cu = 0.04, -0.015, 0.02
x = 0.5 + np.array([-4, -2, -1, 0, 1, 2, 4]) * EPS / 2
t = 1.5 * np.concatenate([[0.1, 0.3, 0.49], x, [0.51, 0.7, 2.0]])
y = curve(np.concatenate([t, t]), np.repeat([lv, 2 * lv], len(t)), np.repeat([sl, 2 * sl], len(t)),
          np.repeat([cu, 2 * cu], len(t)))
one, two = y[:len(t)], y[len(t):]
assert np.all(np.abs(one - closed(t, lv, sl, cu)) <= 4 * EPS * (lv + abs(sl) + cu))
seam = one[3:10]
assert np.max(seam) - np.min(seam) <= 4 * EPS * (lv + abs(sl) + cu), seam
assert np.array_equal(two, 2 * one)
print('across x = 1/2:', seam)

lv, sl, cu = 0.04, -0.015, 0.02 x = 0.5 + np.array([-4, -2, -1, 0, 1, 2, 4]) * EPS / 2 t = 1.5 * np.concatenate([[0.1, 0.3, 0.49], x, [0.51, 0.7, 2.0]]) y = curve(np.concatenate([t, t]), np.repeat([lv, 2 * lv], len(t)), np.repeat([sl, 2 * sl], len(t)), np.repeat([cu, 2 * cu], len(t))) one, two = y[:len(t)], y[len(t):] assert np.all(np.abs(one - closed(t, lv, sl, cu)) <= 4 * EPS * (lv + abs(sl) + cu)) seam = one[3:10] assert np.max(seam) - np.min(seam) <= 4 * EPS * (lv + abs(sl) + cu), seam assert np.array_equal(two, 2 * one) print('across x = 1/2:', seam)

across x = 1/2: [0.03180408 0.03180408 0.03180408 0.03180408 0.03180408 0.03180408
 0.03180408]