Case study: fitting a curve by least squares

The Nelson–Siegel curve

Write $y(m)$ for the zero rate, in percent, of maturity $m$ years (case study 1.6.4). The Nelson–Siegel curve builds it from three shapes:

$g_2$ is largest at $u = 1.79$, where it is 0.2984, so the hump peaks at a maturity of about $1.79\tau$. With $\tau = 2$:

Case study 3.6.5 found that yield changes move mostly in three shapes: level, slope and curvature. Nelson–Siegel builds the curve itself from three such shapes: a constant, $g_1$ falling from 1 towards 0, and the hump $g_2$.

The quotes

The ten quotes are invented for illustration. They come from a Nelson–Siegel curve with $\boldsymbol{\beta} = (4.6, -1.8, 1.5)$ and $\tau = 1.8$, plus noise of typical size 3 basis points. The noise is drawn with the seed np.random.default_rng(445), and each quote is rounded to 0.001%.

A basis point (bp) is 0.01 percentage points (case study 1.6.4).

The fit chooses $\boldsymbol{\beta} = (\beta_0, \beta_1, \beta_2)$ and $\tau$ to make the sum of squared errors smallest:

Here $y_i$ is the quote at maturity $m_i$, and $S$ is in squared percentage points. An $S$ of 0.002374 over ten quotes is a root-mean-square miss of 1.54 bp.

With τ fixed: linear least squares

Fix $\tau$, and the curve is a combination of three known columns: a column of ones, $g_1(m_i/\tau)$ and $g_2(m_i/\tau)$. Stack them in a $10 \times 3$ matrix $X$. Then $S = \lVert\mathbf{y} - X\boldsymbol{\beta}\rVert^2$, exactly lesson 3.5.3's least squares, solved with lstsq.

Lesson 4.4.1 showed that this $S$ is a quadratic in $\boldsymbol{\beta}$ with Hessian $2X^\top X$. The three columns are independent, so $2X^\top X$ is positive definite and the normal equations give the one global minimum.

With $\tau = 2$, $\boldsymbol{\beta} = (4.5457, -1.7217, 1.8093)$ and $S = 0.005053$. The residuals, quote minus fit, are in bp:

The 7-year quote is missed by 5.53 bp. The Hessian's eigenvalues are 0.1128, 1.3376 and 24.0792, all positive.

With τ free: a nonlinear problem

$\tau$ sits inside the exponentials, so $S$ is not a quadratic in all four numbers, and $\nabla S = \mathbf{0}$ is not a linear system. This is the problem lessons 4.4.2 and 4.4.3 were built for.

The profile: one τ at a time

Each $\tau$ has its own best $\boldsymbol{\beta}$, from lstsq. Call the sum of squares it leaves $S^*(\tau)$, the profile. The lowest point of the profile is the lowest point of $S$ over all four numbers.

The profile has two valleys: a local minimum near $\tau = 0.32$ and the global one near $\tau = 1.63$, with a hump near $\tau = 0.69$ between them. These are lesson 4.4.1's local and global minima, met in a real problem.

A grid search tries $\tau = 0.1, 0.11, \ldots, 10$, a linear fit at each, and keeps the best. It finds $\tau = 1.63$, with $\boldsymbol{\beta} = (4.5868, -1.8404, 1.5298)$ and $S = 0.002374$. It cannot be fooled by the second valley, so many desks fit this way; the exercise asks you to write it.

All four numbers at once

To polish the grid's answer, treat all four numbers as one point $\mathbf{p} = (\beta_0, \beta_1, \beta_2, \tau)$ and minimise $S(\mathbf{p})$. Its gradient is $2\sum_i (y(m_i) - y_i)\,\nabla y(m_i)$.

The partial derivatives of $y(m)$ in the three $\beta$s are 1, $g_1$ and $g_2$. The one in $\tau$ comes from the chain rule (lesson 2.3.3), as in the code.

The cells take the Hessian by central differences of the gradient (lesson 4.1.2), with a gap of $10^{-5}$. Two Newton steps from the grid's answer reach the minimum:

There the Hessian has eigenvalues 0.0196, 0.1796, 1.5573 and 24.0767. All are positive, so this is a strict local minimum (lesson 4.4.1), and the profile shows it is the global one. Their ratio, the condition number, is about 1229: a long narrow valley (lesson 4.4.2).

Gradient descent must keep $\eta$ below $2/24.0767 = 0.0831$. From the $\tau = 1$ fit with $\eta = 0.08$, it is still 0.3839 from the minimum after 1,000 steps. It needs 5,000 steps to come within 0.0009.

Plain Newton from the same $\tau = 1$ fit goes badly. There the Hessian has an eigenvalue of $-0.0028$, so it is not positive definite (lesson 4.4.3). The step lands on $\tau = -0.3158$, a negative time scale with no meaning.

Started from the grid's answer, where the Hessian is positive definite, Newton takes two steps.

scipy.optimize.minimize, and the other valley

minimize(sse, start) (lesson 4.4.3) runs BFGS on all four numbers. From the start $(4, -1, 1, 1)$ it takes 20 iterations and returns the best fit.

From $(5, 0, 0, 1)$ it takes 36 iterations, reports success, and returns a different curve. The counts can differ slightly with the SciPy version.

The other fit is the profile's local minimum near $\tau = 0.32$. Its $S$ is 3.12 times larger, and its $\beta_2$ is negative: a dip near the short end instead of a hump. The two curves are close across the quotes but part company before the first one, at 6 months.

Parameters against the curve. The quotes came from $\tau = 1.8$ and $\boldsymbol{\beta} = (4.6, -1.8, 1.5)$, yet the best fit has $\tau = 1.6271$. Its curve still stays within 3.48 bp of the true curve from 3 months to 30 years.

Noise moves the parameters easily along the flat directions of $S$, those with small Hessian eigenvalues, while the curve barely moves. So read the curve, not the individual numbers.

On the desk

import numpy as np

MATURITIES = np.array([0.5, 1, 2, 3, 5, 7, 10, 15, 20, 30], dtype=float)
QUOTES = np.array([3.2, 3.517, 3.953, 4.213, 4.438, 4.46, 4.538, 4.559, 4.578, 4.556])
GRID = np.round(np.arange(0.1, 10.0001, 0.01), 2)


def _expect_value_error(call, why):
    try:
        call()
    except ValueError:
        return
    raise AssertionError(f"expected ValueError: {why}")


def test_curve_values():
    """ns_curve with beta = (4, −1, 1) and tau = 2 gives 3.3935 at 1 year and 3.9933 at 10 years"""
    got = ns_curve([1.0, 10.0], (4.0, -1.0, 1.0), 2.0)
    expect_close(got, np.array([3.393469340287367, 3.9932620530009144]), what="ns_curve([1, 10], (4, -1, 1), 2)")


def test_curve_long_end():
    """Far out, the curve tends to beta0: 4.0 at a million years"""
    expect_close(ns_curve([1e6], (4.0, -1.0, 1.0), 2.0), np.array([4.0]), what="ns_curve at 1,000,000 years", rel=1e-5)


def test_case_study_fit():
    """The case study's quotes on the grid 0.1, 0.11, ..., 10 give tau = 1.63 and S = 0.002374"""
    beta, tau, sse = fit_nelson_siegel(MATURITIES, QUOTES, GRID)
    expect_close(tau, 1.63, what="tau")
    expect_close(beta, np.array([4.586793133757306, -1.8404383921777805, 1.529770785577442]), what="beta", rel=1e-6)
    expect_close(sse, 0.0023738896016146332, what="the sum of squared errors", rel=1e-6)


def test_return_types():
    """It returns beta as a NumPy array of 3 numbers, then tau and the sum of squares as numbers"""
    beta, tau, sse = fit_nelson_siegel(MATURITIES, QUOTES, [1.0, 2.0])
    assert isinstance(beta, np.ndarray) and beta.shape == (3,), f"beta should be an array of 3 numbers, got {beta!r}"
    assert np.isscalar(tau) and np.isscalar(sse), f"tau and sse should be numbers, got {tau!r} and {sse!r}"


def test_exact_curve():
    """Quotes made from beta = (5, −2, 1) and tau = 3, with 3 on the grid, are fitted exactly"""
    m = np.array([0.25, 1, 2, 4, 6, 8, 12, 25])
    u = m / 3.0
    g1 = (1 - np.exp(-u)) / u
    y = 5.0 - 2.0 * g1 + 1.0 * (g1 - np.exp(-u))
    beta, tau, sse = fit_nelson_siegel(m, y, [1.0, 2.0, 3.0, 4.0])
    expect_close(tau, 3.0, what="tau")
    expect_close(beta, np.array([5.0, -2.0, 1.0]), what="beta")
    expect_close(sse, 0.0, what="the sum of squared errors", abs_tol=1e-12)


def test_grid_decides_the_valley():
    """A grid that stops below 0.6 finds the other valley: tau = 0.32 and S = 0.007409"""
    beta, tau, sse = fit_nelson_siegel(MATURITIES, QUOTES, np.round(np.arange(0.1, 0.6, 0.01), 2))
    expect_close(tau, 0.32, what="tau")
    expect_close(sse, 0.007409477119638619, what="the sum of squared errors", rel=1e-6)


def test_lists_work():
    """Plain Python lists work as inputs"""
    beta, tau, sse = fit_nelson_siegel(MATURITIES.tolist(), QUOTES.tolist(), [1.63])
    expect_close(sse, 0.0023738896016146332, what="the sum of squared errors at tau = 1.63", rel=1e-6)


def test_lengths_differ():
    """Ten maturities and nine yields: ValueError"""
    _expect_value_error(lambda: fit_nelson_siegel(MATURITIES, QUOTES[:9], GRID), "the inputs have different lengths")


def test_too_few_quotes():
    """Three quotes cannot fix four parameters: ValueError"""
    _expect_value_error(lambda: fit_nelson_siegel(MATURITIES[:3], QUOTES[:3], GRID), "3 quotes for 4 parameters")


def test_missing_quote():
    """A missing quote (nan) is refused: ValueError"""
    bad = QUOTES.copy()
    bad[4] = np.nan
    _expect_value_error(lambda: fit_nelson_siegel(MATURITIES, bad, GRID), "the 5-year quote is nan")


def test_bad_maturity_or_tau():
    """A zero maturity, a negative tau or an empty grid: ValueError"""
    zero = MATURITIES.copy()
    zero[0] = 0.0
    _expect_value_error(lambda: fit_nelson_siegel(zero, QUOTES, GRID), "a maturity of 0 divides by zero in g1")
    _expect_value_error(lambda: fit_nelson_siegel(MATURITIES, QUOTES, [1.0, -2.0]), "tau = -2 is not positive")
    _expect_value_error(lambda: fit_nelson_siegel(MATURITIES, QUOTES, []), "the grid of taus is empty")