Calculate bond yield to maturity (YTM) in C++

Question

I want to calculate bond's yield to maturity given price using either the bisection or the secant method. I know there's C++ recipes online for this, but I can't figure out what's wrong with my own code. Both methods create an infinite loop. When I tried to step through variables, the yield to maturity iterative solutions keep increasing to 80% and higher.

#include <iostream>
#include <math.h>

using namespace std;

class bond {
private:
    double principal;
    double coupon;
    double timeToMaturity;

public:
    bond(double principal, double coupon, double timeToMaturity) {
        this->principal = principal;
        this->coupon = coupon;
        this->timeToMaturity = timeToMaturity;
    }

    double getprice(double YTM);
    double getytm(double price);
    double getytmsec(double price);
};

bool isZero(double x)
{
    if (fabs(x) < 1e-10)
        return true;
    return false;
}

double bond::getprice(double YTM) {
    double pvPrincipal;
    double pvCoupon;
    double factor = 0;

    pvPrincipal = this->principal / pow(1 + YTM, this->timeToMaturity);

    for (int i = 0; (this->timeToMaturity - i) > 0; i++) {
        double denom = pow(1 + YTM, this->timeToMaturity - i);
        factor += 1 / denom;
    }

    pvCoupon = this->coupon * factor;

    return pvPrincipal + pvCoupon;
}

// Bisection method
double bond::getytm(double price) {
    double low = 0;
    double high = 1;

    double f0 = getprice(low);
    double f2 = 1;
    double x2 = 0;

    while (!isZero(f2)) {
        x2 = (low + high) / 2;
        f2 = getprice(x2) - price;
        if (f2 < 0) {
            low = x2;
        }
        else {
            high = x2;
        }
    }

    return x2;
}

// Secant method
double bond::getytmsec(double price) {
    double x1 = price;
    double x2 = price + 0.25;
    double f1 = this->getprice(x1);
    double f2 = this->getprice(x2);
    for (; !isZero(f1 - price); x2 = x1, x1 = price) {
        f1 = getprice(x1);
        f2 = getprice(x2);
        price = x1 - f1 * (x1 - x2) / (f1 - f2);
    }
    return price;
}

int main() {
    bond myBond = { 1000, 25, 6 };
    cout << "YTM is " << myBond.getytm(950) << endl;
    cout << "YTM is " << myBond.getytmsec(950) << endl;

    return 0;
}

Answer 1

As suggested, a good way to debug this is to step through the calculations. Alternatively, you can print the relevant values at each iteration.

The problem is: find a zero to the function f(x) = getprice(x) - price.

The bisection method in general is: start with an interval [low, high] where f(low) and f(high) are of different signs (one non-positive, one non-negative). That means it contains a zero. Then choose either the left or right sub-interval based on the function value at the midpoint to maintain that property.

In this case, the function is monotonic and non-increasing so we know that f(low) must be the larger (non-negative) number and f(high) must be the smaller (non-positive) number. Therefore, we must choose the left sub-interval if f(midpoint) is negative, and choose the right sub-interval if f(midpoint) is positive. But the code does the opposite, choosing the right sub-interval if f(midpoint) is negative:

    x2 = (low + high) / 2;
    f2 = getprice(x2) - price;
    if (f2 < 0) {
        low = x2;
    }

So you choose smaller and smaller right sub-intervals with eventually [low, high] = [1, 1] and it is an infinite loop. Replace f2 < 0 with f2 > 0.

The secant method generally involves taking two zero "estimates" x_k and x_{k-1} and uses a recurrence to find a better "estimate" x_{k+1}. The recurrence essentially uses the line between (x_{k-1}, f(x_{k-1}) and (x_k, f(x_k)) and looks where this line crosses zero.

The code provided has multiple problems. First, in the important step:

    price = x1 - f1 * (x1 - x2) / (f1 - f2);

where x1 and x2 are the current and previous estimates and f1 is getprice(x1) and f2 is getprice(x2). Importantly, note that f1 is not f(x1) where f is the function whose zero we want. This is not the secant formula. The first part of the second term should be the value of the function at x1, i.e. f1 - price, and not f1:

    ... = x1 - (f1 - price) * (x1 - x2) / (f1 - f2);

Secondly, you are assigning this to price and thus losing the actual value of price which you do need at each iteration.

Thirdly, the initial guesses for yield are price and price + 0.25. These are so far away from the actual value that it becomes a problem (the zero is a yield, between 0 and 1). Try 0 and 1.

A lot of this could be avoided by not intermixing many concerns. You could factor out the logic of finding a zero of a function from the actual identity of the function. For example, one step of bisection is:

template<typename Function>
constexpr auto bisection_step(double low, double high, Function f) -> std::pair<double, double>
{
    assert(std::isfinite(high));
    assert(std::isfinite(low));
    assert(low < high);
    assert(f(low) * f(high) <= 0.);

    auto mid = midpoint(low, high);
    if (f(low) * f(mid) <= 0.)
        return {low, mid};
    else
        return {mid, high};
}

That lets you specify assumptions in the form of assertions or checks that throw exceptions or return error codes. That also makes the logic clearer, so one is less likely to select the wrong sub-interval. And even if one did, the assertion would fire.