Modified duration¶

A bond's modified duration, the relative fall of its price per unit of yield, D = -(1/P) dP/dy, by the central difference of its prices a bump either side of its yield, D = (P(y - h) - P(y + h)) / (2 h P); for a bullet bond, Macaulay's duration over (1 + y/f).

The prices are a game's, lusine.finance.bond_price's at the three yields, each per unit of face, so in no currency. Source: Fabozzi, Bond Markets, Analysis and Strategies, ch. 4. Checked against a reference in numpy, which has its own tests.

Source

Fabozzi, Bond Markets, Analysis and Strategies, ch. 4.

Proof

curve, limit.

duration.trait

namespace lusine.finance;

float<year> duration(float price, float price_down, float price_up, float<year^-1> bump)
{
  return (price_down - price_up) / (2.0 * bump * price);
}

Uses, without proving it here: finance.bond_price.

A bond's modified duration, the relative fall of its price per unit of yield, by the central difference of its prices a bump either side of its yield. The prices are the library's lusine.finance.bond_price, run first at the three yields, as a game would run it before this kernel.

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/duration.trait'))

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

// lusine.finance.duration: a bond's modified duration, the relative fall of its price per unit of yield,
// D = -(1/P) dP/dy, by the central difference of its prices a bump either side of its yield,
// D = (P(y - h) - P(y + h)) / (2 h P); for a bullet bond, Macaulay's duration over (1 + y/f). The prices are a game's,
// lusine.finance.bond_price's at the three yields, each per unit of face, so in no currency.
// Source: Fabozzi, Bond Markets, Analysis and Strategies, ch. 4.
// Checked against a reference in numpy, which has its own tests.

namespace lusine.finance;

float<year> duration(float price, float price_down, float price_up, float<year^-1> bump)
{
  return (price_down - price_up) / (2.0 * bump * price);
}

The program¶

Two programs, each a lattice of items, one a case: the first runs bond_price at each bond's yield and a bump either side, the second this kernel on those prices. 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 the duration, in years, is read back over YEAR. macaulay is the independent account, the cash flows' present-value-weighted time over (1 + y/f), in double precision. The bump is a basis point. tolerance is binary32's: each price within the bond page's bound of its value, the difference of two divided by 2 h, and the central difference's own error, below a part in 10^5.

In [3]:

print(trait_run.law_program('finance/bond_price.trait', 'bond_price'))
print(trait_run.law_program('finance/duration.trait', 'duration'))

def prices(coupon, years, yld, bump, freq=2.0):
    """Bonds' prices per 100 at their yields and a bump either side, by the library's bond_price: one lattice."""
    coupon, years, yld = np.broadcast_arrays(*(np.asarray(a, dtype=float) for a in (coupon, years, yld)))
    at = np.concatenate([yld, yld - bump, yld + bump])
    p = trait_run.law('finance/bond_price.trait', 'bond_price', face=100.0, coupon=np.tile(coupon, 3) / YEAR,
                      years=np.tile(years, 3) * YEAR, yld=at / YEAR, freq=freq / YEAR)
    return np.split(p, 3)


def duration(price, price_down, price_up, bump):
    """The kernel's modified duration, in years."""
    return trait_run.law('finance/duration.trait', 'duration', price=price, price_down=price_down,
                         price_up=price_up, bump=np.asarray(bump) / YEAR) / YEAR


def macaulay(coupon, years, yld, freq=2.0):
    """The independent account: Macaulay's duration, each cash flow's time weighted by its present value, over
    (1 + y/f), in double precision."""
    out = []
    for c, n, y in np.broadcast(coupon, years, yld):
        k = np.arange(1, int(round(n * freq)) + 1)
        flows = np.full(len(k), 100.0 * c / freq)
        flows[-1] += 100.0
        pv = flows * (1 + y / freq) ** -k.astype(float)
        out.append(np.sum(k / freq * pv) / np.sum(pv) / (1 + y / freq))
    return np.array(out)


def tolerance(coupon, years, yld, bump, want, freq=2.0):
    """In years: each price's binary32 bound, the bond page's (8 ulps, and the rounding of 1 + y carried n periods
    and through the annuity), twice over 2 h; and the central difference's own error, below a part in 10^5."""
    n, y = np.round(np.asarray(years) * freq), np.asarray(yld) / freq
    c, v = 100.0 * np.asarray(coupon) / freq, (1 + y) ** -n
    price = c * (1 - v) / y + 100.0 * v
    return 2 * EPS * (8 + n / 2 + c * n * v / (2 * y * price)) / (2 * bump) + 1e-5 * want

print(trait_run.law_program('finance/bond_price.trait', 'bond_price')) print(trait_run.law_program('finance/duration.trait', 'duration')) def prices(coupon, years, yld, bump, freq=2.0): """Bonds' prices per 100 at their yields and a bump either side, by the library's bond_price: one lattice.""" coupon, years, yld = np.broadcast_arrays(*(np.asarray(a, dtype=float) for a in (coupon, years, yld))) at = np.concatenate([yld, yld - bump, yld + bump]) p = trait_run.law('finance/bond_price.trait', 'bond_price', face=100.0, coupon=np.tile(coupon, 3) / YEAR, years=np.tile(years, 3) * YEAR, yld=at / YEAR, freq=freq / YEAR) return np.split(p, 3) def duration(price, price_down, price_up, bump): """The kernel's modified duration, in years.""" return trait_run.law('finance/duration.trait', 'duration', price=price, price_down=price_down, price_up=price_up, bump=np.asarray(bump) / YEAR) / YEAR def macaulay(coupon, years, yld, freq=2.0): """The independent account: Macaulay's duration, each cash flow's time weighted by its present value, over (1 + y/f), in double precision.""" out = [] for c, n, y in np.broadcast(coupon, years, yld): k = np.arange(1, int(round(n * freq)) + 1) flows = np.full(len(k), 100.0 * c / freq) flows[-1] += 100.0 pv = flows * (1 + y / freq) ** -k.astype(float) out.append(np.sum(k / freq * pv) / np.sum(pv) / (1 + y / freq)) return np.array(out) def tolerance(coupon, years, yld, bump, want, freq=2.0): """In years: each price's binary32 bound, the bond page's (8 ulps, and the rounding of 1 + y carried n periods and through the annuity), twice over 2 h; and the central difference's own error, below a part in 10^5.""" n, y = np.round(np.asarray(years) * freq), np.asarray(yld) / freq c, v = 100.0 * np.asarray(coupon) / freq, (1 + y) ** -n price = c * (1 - v) / y + 100.0 * v return 2 * EPS * (8 + n / 2 + c * n * v / (2 * y * price)) / (2 * bump) + 1e-5 * want

namespace probe;

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

struct item {
  host float face;
  host float<year^-1> coupon;
  host float<year> years;
  host float<year^-1> yld;
  host float<year^-1> freq;
  float result_0;
};
lattice<item, 1> items(spacing = 1 m, boundary = wall);

void main() { items.result_0 = lusine.finance.bond_price(items.face, items.coupon, items.years, items.yld, items.freq); }

namespace probe;

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

struct item {
  host float price;
  host float price_down;
  host float price_up;
  host float<year^-1> bump;
  float<year> result_0;
};
lattice<item, 1> items(spacing = 1 m, boundary = wall);

void main() { items.result_0 = lusine.finance.duration(items.price, items.price_down, items.price_up, items.bump); }

Curve: duration grows with the tenor, less than the tenor¶

A 5% bond at 4%, from one year to thirty, against Macaulay's over (1 + y/f): rising, and below the tenor.

In [4]:

years, bump = np.arange(1.0, 31.0), 1e-4
d = duration(*prices(0.05, years, 0.04, bump), bump)
want = macaulay(0.05, years, 0.04)
tol = tolerance(0.05, years, 0.04, bump, want)
assert np.all(np.abs(d - want) <= tol), np.max(np.abs(d - want) / tol)
assert np.all(np.diff(d) > 0) and np.all(d < years)
fig, ax = plt.subplots(figsize=(7, 3))
ax.plot(years, d, 'o', ms=4, label='duration')
ax.plot(years, want, lw=0.8, label="Macaulay's / (1 + y/f)")
ax.plot(years, years, '--', lw=0.8, label='tenor')
ax.set_xlabel('tenor (years)')
ax.set_ylabel('years')
ax.legend()
plt.show()

years, bump = np.arange(1.0, 31.0), 1e-4 d = duration(*prices(0.05, years, 0.04, bump), bump) want = macaulay(0.05, years, 0.04) tol = tolerance(0.05, years, 0.04, bump, want) assert np.all(np.abs(d - want) <= tol), np.max(np.abs(d - want) / tol) assert np.all(np.diff(d) > 0) and np.all(d < years) fig, ax = plt.subplots(figsize=(7, 3)) ax.plot(years, d, 'o', ms=4, label='duration') ax.plot(years, want, lw=0.8, label="Macaulay's / (1 + y/f)") ax.plot(years, years, '--', lw=0.8, label='tenor') ax.set_xlabel('tenor (years)') ax.set_ylabel('years') ax.legend() plt.show()

Limit: a zero-coupon bond's is its tenor over (1 + y/f)¶

With no coupons there is one cash flow, and its Macaulay duration is the tenor itself. A price that does not move has no duration.

In [5]:

years, bump = np.array([1.0, 5.0, 10.0, 30.0]), 1e-4
d = duration(*prices(0.0, years, 0.04, bump), bump)
want = years / (1 + 0.04 / 2)
assert np.all(np.abs(d - want) <= tolerance(0.0, years, 0.04, bump, want)), d - want
flat = duration(np.array([97.0]), np.array([97.0]), np.array([97.0]), bump)
assert flat[0] == 0.0
print('zero coupon:', d, 'against', want)

years, bump = np.array([1.0, 5.0, 10.0, 30.0]), 1e-4 d = duration(*prices(0.0, years, 0.04, bump), bump) want = years / (1 + 0.04 / 2) assert np.all(np.abs(d - want) <= tolerance(0.0, years, 0.04, bump, want)), d - want flat = duration(np.array([97.0]), np.array([97.0]), np.array([97.0]), bump) assert flat[0] == 0.0 print('zero coupon:', d, 'against', want)

zero coupon: [ 0.98069315  4.90259437  9.80526377 29.41665323] against [ 0.98039216  4.90196078  9.80392157 29.41176471]

Invariant: a ratio, in no currency¶

The duration is the prices' relative fall, so prices per 200 of face, every one doubled, give the same duration bit for bit.

In [6]:

bump = 1e-4
p, down, up = prices(0.05, np.array([2.0, 10.0, 30.0]), 0.04, bump)
d = duration(np.concatenate([p, 2 * p]), np.concatenate([down, 2 * down]), np.concatenate([up, 2 * up]), bump)
assert np.array_equal(d[:3], d[3:]), d
print('per 100 and per 200:', d[:3])

bump = 1e-4 p, down, up = prices(0.05, np.array([2.0, 10.0, 30.0]), 0.04, bump) d = duration(np.concatenate([p, 2 * p]), np.concatenate([down, 2 * down]), np.concatenate([up, 2 * up]), bump) assert np.array_equal(d[:3], d[3:]), d print('per 100 and per 200:', d[:3])

per 100 and per 200: [ 1.89117956  7.92203514 16.59897483]