C++ bond mathematics calculator

Question

I have written a program in C++ that processes financial bond mathematics, making extensive use of the valarray class, so that mathematical functions and operations can be applied vectorially on whole arrays in just one expression, instead of having to resort to a loop over all elements.

This is similar to how numeric data is processed in languages more focused on scientific computation like Fortran or Matlab.

There are many function templates in my code that perform one-line operations as described above, which is ideal if the template parameter is instantiated as a scalar, say an int or double, or as a valarray<int or double>. I guess this is less ideal instead if the template parameter is taken to be an array or a vector. This perhaps might constitute a possible objection that can be levelled to the logic of my program, but I don't think it's practical to introduce loops just so one can accept vectors as data.

The mathematics behind the program revolve around the equation for the bond value \$B\$:

$$ B = \sum_{i=1}^{n-1} \frac{C}{2} P \cdot d(t_i) + \left( 1 + \frac{C}{2} \right) P \cdot d(t_n) \qquad\qquad\qquad\qquad\qquad\qquad (1) $$

where \$t_i\$ are the \$i = 1, 2 \dots n\$ coupon dates, \$C\$ is the bond coupon rate, \$P\$ is the bond principal (initial loan to be returned at maturity \$t_n\$), \$d(t_i)\$ are the discount factors for each coupon date.

The formula to obtain \$d(t_i)\$ varies based on the type of interest rate \$r(t_i)\$ considered:

• Zero-rate: $$d(t_i) = e^{- r(t_i) \cdot t_i} \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad (2)$$ where \$r(t_i)\$ is given either tabular or as a continuous function of \$t\$.
• Instantaneous rate: $$d(t_i) = e^{- \int_{0}^{t_i} r(\tau) d\tau} \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad (3)$$ where \$r(t)\$ is a given continuous function.

Other bond-related quantities the program computes are:

• the yield \$y\$, which is the constant interest rate that substituted in (1) attains a known \$B\$. That is given \$B\$, solve (1) implicitly to extract \$y\$ with a numerical method (I hardcoded Newton-Raphson).
• The par yield coupon rate is the value of \$C\$ that makes \$B = P\$ in (1).
• The bond duration \$D\$ is the opposite of the relative change of \$B\$ coming after a change of yield \$y\$: $$D = - \frac{1}{B} \frac{\partial B}{\partial y}$$
• The bond convexity \$V\$ is the opposite of the relative change of duration with respect to yield: $$V = \frac{1}{B} \frac{\partial^2 B}{\partial y^2}$$

The program below creates a Bond class hierarchy with four derived classes:

• KnownRateBond takes zero-interest rate input in tabular form;
• ZeroRateBond takes zero-interest rate input as a continuous function of time;
• InstantRateBond takes instantaneous interest rate input, so it has to evaluate the integral in (3) with NewtonCotesFormulas classes of quadratures that I had written here;
• BondGivenYield takes the yield value as input instead of the interest rate to compute the other quantities.

The code is below, please don't spare any critique or suggestion for possible enhancements:

main.cpp:

/** Computes the important quantities for a bond from input
*   All quantities are acquired as class data members,
*   useful to print all on file or on excel
*   TODO: - add parser to type interest time-law directly on input
*         - add xlw interface
*         - fix virtual computeDiscountFactors() function not overridden by InstantRateBond class
*/
#include <iostream>
#include <initializer_list>
#include <algorithm>
#include <iterator>
#include <cmath>
#include <valarray>
#include <functional>
#include "Bond.h"

using std::initializer_list;    using std::valarray;
using std::cout;                using std::cin;
using std::function;            using std::copy;
using std::begin;               using std::unique_ptr;

// helps with input
void flushCin() {
    cin.ignore();
    cin.clear();
}

// processes inputs
void inputFromUser(valarray<double>& times,
                   double& coupon,
                   valarray<double>& inp_known_rates,
                   function<valarray<double> (valarray<double>)>& inp_rate,
                   double& inp_yield,
                   InterestRates& type,
                   Formula& quadrature,
                   double& principal,
                   double& tol)
{
    cout << "Please input coupon dates as fraction (decimal number) of the year,"
         << " and press Enter key when done: \n";
    constexpr size_t MAX_SIZE = 1000;
    double input_times[MAX_SIZE];
    int counter;
    while (cin.peek() != '\n')  // while next input is not Enter key
        cin >> input_times[counter++];
    flushCin();
    times.resize(counter);
    copy(input_times, input_times+counter, begin(times));
    // or briefly, while testing
    //times = {0.1667, 0.6667, 1.1667, 1.6667};
    //times = {0.5, 1., 1.5, 2};
    inp_known_rates.resize(times.size());

    cout << "\nPlease enter the coupon rate in fraction of year (decimal): \n";
    cin >> coupon;

    cout << "\nPlease select (1, 2, 3 or 4) how the interest rate law above should be computed: \n"
         << "Known rates at given dates ---> (1)\n"
         << "Zero rate --------------------> (2)\n"
         << "Instantaneous rate -----------> (3)\n"
         << "Yield rate -------------------> (4)\n";
    int input_type;
    do {
        cin >> input_type;
        switch (input_type) {
            case 1:
               type = InterestRates::known_rate;
               break;
            case 2:
               type = InterestRates::zero_rate;
               break;
            case 3:
               type = InterestRates::instantaneous;
               break;
            case 4:
                type = InterestRates::given_yield;
                break;
            default:
                cout << "Select 1, 2, 3 or 4 only\n";
                flushCin();
                break;
        }
    } while (input_type < 1 || input_type > 4);
    flushCin();

    cout << "\nPlease digit a value for the principal, or press Enter to leave it at 100 : \n";
    while (cin.peek() != '\n')
        cin >> principal;
    flushCin();

    if (type == InterestRates::known_rate) {
//        inp_known_rates = {0.05, 0.0525, 0.0535, 0.055};  // debug
        cout << "Please input known interest rate at coupon dates, separated by space,"
             << " and press Enter key when done: \n";
        while (cin.peek() != '\n') {  // while next input is not Enter key
            for (auto& rate : inp_known_rates)
                cin >> rate;
        }
        flushCin();
    } else if (type == InterestRates::instantaneous) {
        cout << "Choose the numerical integration method for the interest rate law :\n"
             << "Midpoint numeric quadrature -----> (1)\n"
             << "Trapezoidal numeric quadrature --> (2)\n"
             << "Simpson's numeric quadrature ----> (3)\n\n";
        int input_quadrature;
        do {
            cin >> input_quadrature;
            switch (input_quadrature) {
                case 1:
                    quadrature = Formula::Midpoint;
                    break;
                case 2:
                    quadrature = Formula::Trapezoidal;
                    break;
                case 3:
                    quadrature = Formula::Simpsons;
                    break;
                default:
                    cout << "Select 1, 2 or 3 only\n";
                    flushCin();
                    break;
            }
        } while (input_quadrature < 1 || input_quadrature > 3);
        flushCin();

        cout << "\nPlease enter the tolerance value sufficient for convergence of integration above,";
        cout << " or press Enter to leave it at 1E-6 : \n";
        while (cin.peek() != '\n')
            cin >> tol;
        flushCin();
    } else if (type == InterestRates::given_yield) {
        cout << "Please input the bond yield and press Enter key when done: \n";
        while (cin.peek() != '\n') {  // while next input is not Enter key
            cin >> inp_yield;
        }
    }
}

// functor for interest rate law if any
// would be nice to input this directly through a parser
template<typename T>
function<T (T)> interestRateLaw() {
    return [](const T t)->T {
               //return 0.0525 + 1 / (100 * (1 + exp(- t * t)));
               return 0.0525 + log(1 + 2.*t) / 200;  // another law
    };
}


int main()
{
    // declare variables
    valarray<double> times, inp_known_rates;
    double coupon, principal{100}, tol{1E-6}, inp_yield{0};
    InterestRates bond_type;
    Formula quadrature;
    // acquire input
    function<valarray<double> (valarray<double>)> inp_rate = interestRateLaw<valarray<double> >();
    inputFromUser(times, coupon, inp_known_rates, inp_rate, inp_yield, bond_type, quadrature, principal, tol);

    // calculates bond data and print it on output
    unique_ptr<Bond> my_bond = makeBond(bond_type, times, coupon, inp_known_rates,
                                        inp_rate, inp_yield, quadrature, principal, tol);

    my_bond->printData(cout, bond_type);

    return 0;
}

bond.h:

#ifndef BOND_H_INCLUDED
#define BOND_H_INCLUDED

#include <iostream>
#include <memory>
#include <functional>
#include <valarray>
#include "NewtonCotesFormulas.h"

// bond types
enum class InterestRates {
    known_rate = 1,
    zero_rate,
    instantaneous,
    given_yield
};

// BASE CLASS
class Bond {
public:
    Bond(std::valarray<double>, const double);
    virtual ~Bond() = default;

    void setConstants(size_t&, double&);

    virtual void computeDiscountFactors(){} // will be overridden by all derived classes except one
    void computeValues(const double, const double = 100, const double = 1E-6);

    std::valarray<double> getCashTimes() const;
    double getCouponRate() const;
    std::valarray<double> getDiscountFactors() const;
    double getValue() const;
    double getDuration() const;
    double getConvexity() const;

    virtual double getGivenYield() const {return 0;}
    double getYield() const;
    double getParYieldCoupon() const;
    double NewtonMethodForYield(const double, const double,
                                std::function<double (double)>&, std::function<double (double)>&);

    void printData(std::ostream&, const InterestRates&);

protected:
    // acquired on input
    std::valarray<double> cash_times;
    double coupon_rate;
    // computed after initialization
    std::valarray<double> discount_factors;
    double bond_value;
    double duration;
    double convexity;
    double yield;
    double par_yield_coupon;
};

// DERIVED CLASSES
class KnownRateBond : public Bond {
public:
    KnownRateBond(std::valarray<double>, const double, std::valarray<double>, double = 100);
    ~KnownRateBond() = default;
    void computeDiscountFactors() override;
private:
    std::valarray<double> known_rates;
};

template<typename T>
class ZeroRateBond : public Bond {
public:
    ZeroRateBond(std::valarray<double>, const double, std::function<T (T)>&, double = 100);
    ~ZeroRateBond() = default;
    void computeDiscountFactors() override;
private:
    std::function<T (T)> interest_rate;
};

template<typename T>
class InstantRateBond : public Bond {
public:
    InstantRateBond(std::valarray<double>, const double, std::function<T (T)>&,
                    const Formula&, double, const double=1E-6);
    ~InstantRateBond() = default;
    void computeDiscountFactors(const Formula&, const double); // alas doesn't override virtual function
private:
    std::function<T (T)> interest_rate;
};

class BondGivenYield : public Bond {
public:
    BondGivenYield(std::valarray<double>, const double, double, double);
    ~BondGivenYield() = default;
    void computeDiscountFactors() override;
    double getGivenYield() const override;
private:
    double given_yield;
};

// BOND FACTORY FUNCTION
template<typename T>
std::unique_ptr<Bond> makeBond(const InterestRates&,
                               std::valarray<double>,
                               const double,
                               std::valarray<double>,
                               std::function<T (T)>&,
                               double,
                               const Formula&,
                               const double = 100,
                               double = 1E-6);

#include "Bond.cpp"
#endif // BOND_H_INCLUDED

bond.cpp:

#include <iostream>
#include <iomanip>
#include <numeric>
#include <cmath>
#include <iterator>
#include <memory>
#include <algorithm>
#include <stdexcept>
#include "Bond.h"

using std::valarray;    using std::function;
using std::begin;       using std::end;
using std::accumulate;  using std::exp;
using std::copy;        using std::setprecision;
using std::cout;        using std::move;
using std::ostream;     using std::ostream_iterator;
using std::make_unique; using std::runtime_error;
using std::unique_ptr;  using std::fixed;
using std::scientific;

// BASE CLASS
Bond::Bond(valarray<double> inp_times, const double inp_coupon)
    : cash_times{inp_times}, coupon_rate{inp_coupon}, discount_factors(cash_times.size()),
    bond_value{0}, duration{0}, convexity{0}, yield{0}, par_yield_coupon{0}
{}

// sets auxiliary variables that are not included in the class as data
void Bond::setConstants(size_t& index_maturity, double& coupon_frequency) {
    if (cash_times.size() == 1) {
        index_maturity = 0;
        coupon_frequency = 1 / cash_times[index_maturity];
    } else {
        index_maturity = cash_times.size() - 1;
        coupon_frequency = 1 / (cash_times[index_maturity] - cash_times[index_maturity - 1]);
    }
}

// compute main bond quantities
void Bond::computeValues(const double coupon_frequency, const double principal, const double tol) {
    // coupon payment
    double coupon_spot_cash_flow = principal * coupon_rate / coupon_frequency;
    // all coupons discounted to present
    valarray<double> cash_flows = coupon_spot_cash_flow * discount_factors;
    // last cashflow at expiry includes return of principal
    *(end(cash_flows)-1) += principal * *(end(discount_factors)-1);
//    copy(begin(cash_flows), end(cash_flows), ostream_iterator<double>(cout, " ")); // debug
//    cout << '\n';

    // main values
    bond_value = accumulate(begin(cash_flows), end(cash_flows), 0.0);
    valarray<double> cash_flows_t = cash_flows * cash_times;
    duration = accumulate(begin(cash_flows_t), end(cash_flows_t), 0.0) / bond_value;
    valarray<double> cash_flows_tt = cash_flows_t * cash_times;
    convexity = accumulate(begin(cash_flows_tt), end(cash_flows_tt), 0.0) / bond_value;

    // auxiliary variables will be used to compute the yield
    valarray<double> exp_cash_times = exp(- cash_times);
    function<double (double)> bond_implicit = [&](double rate_old)
    {
        valarray<double> exp_pow_rate = pow(exp_cash_times, rate_old);
        // *(end(cash_times)-1) below is time of maturity
        return bond_value - principal * exp(- rate_old * *(end(cash_times)-1)) - coupon_spot_cash_flow *
               accumulate(begin(exp_pow_rate), end(exp_pow_rate), 0.0);
    };
    function<double (double)> deriv_bond = [&](double rate_old)
    {
        valarray<double> exp_pow_rate_times_t = cash_times * pow(exp_cash_times, rate_old);
        return principal * *(end(cash_times)-1) * exp(- rate_old * *(end(cash_times)-1)) +
               coupon_spot_cash_flow * accumulate(begin(exp_pow_rate_times_t), end(exp_pow_rate_times_t), 0.0);
    };

    // compute yield and par-yield coupon
    yield = NewtonMethodForYield(tol, 0.1, bond_implicit, deriv_bond);
    par_yield_coupon = (1 - *(end(discount_factors)-1)) * coupon_frequency;
    par_yield_coupon /= accumulate(begin(discount_factors), end(discount_factors), 0.0);
}

// Newton-Raphson method to extract yield
double Bond::NewtonMethodForYield(const double tol, const double rate_init,
                                  function<double (double)>& num, function<double (double)>& denum)
{
    double rate_old = rate_init;
    double rate_new = rate_old - 1;
    while (abs(rate_new - rate_old) > tol) {
        rate_old = rate_new;
        rate_new = rate_old - num(rate_old) / denum(rate_old);
    }

    return rate_new;
}

// getters
valarray<double> Bond::getCashTimes() const {return cash_times;}
double Bond::getCouponRate() const {return coupon_rate;}
valarray<double> Bond::getDiscountFactors() const {return discount_factors;}
double Bond::getValue() const {return bond_value;}
double Bond::getDuration() const {return duration;}
double Bond::getConvexity() const {return convexity;}
double Bond::getYield() const {return yield;}
double Bond::getParYieldCoupon() const {return par_yield_coupon;}

// print all bond data, e.g. on a.txt file
void Bond::printData(ostream& os, const InterestRates& choice) {
    os << "Bond current state :" << setprecision(5) << fixed;

    valarray<double> print_times = getCashTimes();
    os << "\nCoupon dates : ";
    copy(begin(print_times), end(print_times), ostream_iterator<double>(os, " "));
    os << "(per year length)\n";

    valarray<double> print_discounts = getDiscountFactors();
    os << "Discount to present value for each date : ";
    copy(begin(print_discounts), end(print_discounts), ostream_iterator<double>(os, " "));

    os << "\nCoupon rate    : " << setprecision(2) << getCouponRate() * 100 << " % \n";
    os << "Bond VALUE     : " << setprecision(5) << getValue() << '\n';
    os << "Bond duration  : " << getDuration() << '\n';
    os << "Bond convexity : " << getConvexity() << "\n\n";

    if (choice == InterestRates::given_yield)
        os << "Yield given : " << getGivenYield() << '\n';
    os << "Bond yield : " << getYield() << '\n';
    os << "Par yield coupon : " << getParYieldCoupon() << "\n\n";
}


// DERIVATE CLASSES
KnownRateBond::KnownRateBond(valarray<double> inp_times,
                             const double inp_coupon,
                             valarray<double> inp_known_rates,
                             double inp_principal)
                             : Bond{inp_times, inp_coupon}, known_rates{inp_known_rates}
{
    double principal{inp_principal}, coupon_frequency;
//    copy(begin(known_rates), end(known_rates), ostream_iterator<double>(cout, " "));  // debug
//    cout << '\n';
    size_t index_maturity;
    setConstants(index_maturity, coupon_frequency);
    computeDiscountFactors();
    computeValues(coupon_frequency, principal);
}

void KnownRateBond::computeDiscountFactors() {
    discount_factors = exp(- known_rates * cash_times);
}

template<typename T>
ZeroRateBond<T>::ZeroRateBond(valarray<double> inp_times,
                              const double inp_coupon,
                              function<T (T)>& inp_rate,
                              double inp_principal)
                              : Bond(inp_times, inp_coupon), interest_rate(inp_rate)
{
    double principal{inp_principal}, coupon_frequency;
    size_t index_maturity;
    setConstants(index_maturity, coupon_frequency);
    computeDiscountFactors();
    computeValues(coupon_frequency, principal);
}

template<typename T>
void ZeroRateBond<T>::computeDiscountFactors() {
    discount_factors = exp(- interest_rate(cash_times) * cash_times);
}

template<typename T>
InstantRateBond<T>::InstantRateBond(valarray<double> inp_times,
                                    const double inp_coupon,
                                    function<T (T)>& inp_rate,
                                    const Formula& quadrature,
                                    double inp_principal,
                                    const double inp_tol)
                                    : Bond(inp_times, inp_coupon), interest_rate(inp_rate)
{
    double principal{inp_principal}, coupon_frequency, tol{inp_tol};
    size_t index_maturity;
    setConstants(index_maturity, coupon_frequency);
    computeDiscountFactors(quadrature, tol);
    computeValues(coupon_frequency, principal, tol);
}

template<typename T>
void InstantRateBond<T>::computeDiscountFactors(const Formula& quadrature, const double tol) {
    function<T (const T&, const T&, int, function<T (T)>)>
        my_method = createQuadrature<T>(quadrature);
    const size_t array_size = cash_times.size();
    valarray<double> integrated_rate(array_size);
    valarray<double> tol_array(tol, array_size);
    tol_array[array_size-1] *= 1E-2;
    valarray<double> inp_inf(0.0, array_size);
    integrated_rate = convergeToTol(quadrature, inp_inf, cash_times, 4, tol_array, interest_rate, my_method);
    discount_factors = exp(- integrated_rate);
}

BondGivenYield::BondGivenYield(valarray<double> inp_times,
                               const double inp_coupon,
                               double inp_yield,
                               double inp_principal)
                               : Bond{inp_times, inp_coupon}, given_yield{inp_yield}
{
    double principal{inp_principal}, coupon_frequency;
    size_t index_maturity;
    setConstants(index_maturity, coupon_frequency);
    computeDiscountFactors();
    computeValues(coupon_frequency, principal);
}

void BondGivenYield::computeDiscountFactors() {
    discount_factors = exp(- given_yield * cash_times);
}

double BondGivenYield::getGivenYield() const {return given_yield;}

// bond factory function
template<typename T>
unique_ptr<Bond> makeBond(const InterestRates& choice,
                          valarray<double> inp_times,
                          const double inp_coupon,
                          valarray<double> inp_known_rates,
                          function<T (T)>& inp_rate_law,
                          double inp_yield,
                          const Formula& quadrature,
                          const double principal,
                          double tol)
{
    if (choice == InterestRates::known_rate)
        return make_unique<KnownRateBond>(inp_times, inp_coupon, inp_known_rates, principal);
    else if (choice == InterestRates::zero_rate)
        return make_unique<ZeroRateBond<T> >(inp_times, inp_coupon, inp_rate_law, principal);
    else if (choice == InterestRates::instantaneous)
        return make_unique<InstantRateBond<T> >(inp_times, inp_coupon, inp_rate_law,
                                                quadrature, principal, tol);
    else if (choice == InterestRates::given_yield)
        return make_unique<BondGivenYield>(inp_times, inp_coupon, inp_yield, principal);
    else
        throw runtime_error("Pick either choice = 1, 2, 3 or 4");

    return 0;  // never gets here
}

NewtonCotesFormulas.h:

#ifndef NEWTONCOTESFORMULAS_H_INCLUDED
#define NEWTONCOTESFORMULAS_H_INCLUDED

#include <iostream>
#include <functional>

enum class Formula {
    Midpoint = 1,
    Trapezoidal,
    Simpsons
};

std::ostream& operator<< (std::ostream&, const Formula&);

template<typename Quadrature, typename T>
T convergeToTol(const Formula&, const T&, const T&, int, const T&, const std::function<T (T)>&, const Quadrature&);

void printConvergenceValues(std::ostream&, const int, const Formula&, const int, const double, const double, double);

template<typename T> T Midpoint(const T&, const T&, const int, const std::function<T (T)>&);

template<typename T> T Trapezoidal(const T&, const T&, const int, const std::function<T (T)>&);

template<typename T> T Simpsons(const T&, const T&, const int, const std::function<T (T)>&);

template<typename T>
std::function<T (const T&, const T&, const int, const std::function<T (T)>&) > createQuadrature(const Formula&);

#include "NewtonCotesFormulasTemplateImplementation.cpp" // provide definitions for template functions
#endif // NEWTONCOTESFORMULAS_H_INCLUDED

NewtonCotesFormulasTemplateImplementation.cpp:

#include <iostream>
#include <cmath>
#include <functional>
#include <stdexcept>
#include "NewtonCotesFormulas.h"

using std::cout;            using std::abs;
using std::function;        using std::runtime_error;

template<typename Quadrature, typename T>  // also prints convergence table
T convergeToTol(const Formula& choice, const T& inf, const T& sup, int intervals,
                     const T& tol, const function<T (T)>& func, const Quadrature& my_method)
{
    // set up variables
    T old_integral = my_method(inf, sup, intervals, func);
    auto array_size = tol.size();
    printConvergenceValues(cout, 0, choice, intervals, tol[array_size-1],
                           old_integral[array_size-1], old_integral[array_size-1]);
    intervals *= 2;
    T new_integral = my_method(inf, sup, intervals, func);
    T difference = abs(new_integral - old_integral);

    // iterate until convergence
    while((difference > tol).max()) {
        printConvergenceValues(cout, 1, choice, intervals, tol[array_size-1],
                               new_integral[array_size-1], difference[array_size-1]);
        old_integral = new_integral;
        intervals *= 2;
        new_integral = my_method(inf, sup, intervals, func);
        difference = abs(new_integral - old_integral);
    }

    printConvergenceValues(cout, 2, choice, intervals, tol[array_size-1],
                           new_integral[array_size-1], difference[array_size-1]);

    return new_integral;
}

template <typename T>
T Midpoint(const T& inf, const T& sup, const int intervals, const function<T (T) >& func) {
    T interval_width = (sup - inf) / intervals;
    T result(0., inf.size());
    for (int i = 1; i <= intervals; ++i)
        result += func(inf + (i - 0.5) * interval_width);

    return result * interval_width;
}

template <typename T>
T Trapezoidal(const T& inf, const T& sup, const int intervals, const function<T (T)>& func) {
    T interval_width = (sup - inf) / intervals;
    T result = (func(inf) + func(sup)) / 2.;
    for (auto i = 1.; i <= intervals - 1; ++i)
        result += func(inf + i * interval_width);

    return result * interval_width;
}

template <typename T>
T Simpsons(const T& inf, const T& sup, const int intervals, const function<T (T)>& func) {
    T interval_width = (sup - inf) / intervals;
    T result = (func(inf) + func(sup)) / 6;
    for (auto i = 1.; i <= intervals - 1; ++i)
        result += func(inf + i * interval_width) / 3. + 2. * func(inf + (i - 0.5) * interval_width) / 3.;
    result += 2. * func(inf + (intervals - 0.5) * interval_width) / 3.;

    return result * interval_width;
}

// factory function integration methods
template<typename T> function<T (const T&, const T&, const int, const function<T (T)>&) >
createQuadrature(const Formula& choice) {
    if (choice == Formula::Midpoint)
        return Midpoint<T>;
    else if (choice == Formula::Trapezoidal)
        return Trapezoidal<T>;
    else if (choice == Formula::Simpsons)
        return Simpsons<T>;
    else
        throw runtime_error("Pick between choice = 1, 2, or 3");

    return 0;  // never gets here
}

NewtonCotesFormulas.cpp:

#include <iostream>
#include <iomanip>
#include <string>
#include "NewtonCotesFormulas.h"

using std::string;          using std::cout;
using std::setprecision;    using std::fixed;
using std::scientific;      using std::ostream;
using std::to_string;

// this just prints the formula name on cout
ostream& operator<< (ostream& os, const Formula& method) {
    if (method == Formula::Midpoint)
        os << "Midpoint";
    else if (method == Formula::Trapezoidal)
        os << "Trapezoidal";
    else if (method == Formula::Simpsons)
        os << "Simpsons";

    return os;
}

void printConvergenceValues(ostream& os, const int flag, const Formula& choice, const int intervals,
                            const double tol, const double integral, double difference)
{
    string pad_to_size = string(9 - to_string(intervals).size(), ' ');
    switch(flag) {
        case 0:
            os << "Numerical integration, convergence table\n";
            os << choice << " method, with tolerance " << setprecision(1) << tol << "\n";
            os << "intervals" << string(5,' ') << "integral " << string(9,' ') << "tol\n";
            os << setprecision(8) << fixed << pad_to_size << intervals << string(3,' ')
                 << integral << string(3,' ') << scientific << difference << '\n';
            break;
        case 1:
            //cout.setf(ios_base::scientific);
            //cout.precision(8);
            os << pad_to_size << intervals << string(3,' ') << fixed << integral
                 << string(3,' ') << scientific << difference << '\n';
            break;
        case 2:
            os << fixed << pad_to_size << intervals << string(3,' ') << integral
                 << string(3,' ') << scientific << difference << "\n\n";
            break;
    }
}

Results for a semiannual coupon bond with principal \$P = 100\$, coupon rate \$C = 6 \% = 0.06\$, maturity \$t_n = 20\$ months and zero-rate interest law given by \$r(t) = 0.0525 + \frac{\log{(1 + 2 t)}}{200}\$:

Please input coupon dates as fraction (decimal number) of the year, and press En
ter key when done:
0.16667 0.66667 1.16667 1.66667

Please enter the coupon rate in fraction of year (decimal):
0.06

Please select (1, 2, 3 or 4) how the interest rate law above should be computed:

Known rates at given dates ---> (1)
Zero rate --------------------> (2)
Instantaneous rate -----------> (3)
Yield rate -------------------> (4)
2

Please digit a value for the principal, or press Enter to leave it at 100 :

Bond current state :
Coupon dates : 0.16667 0.66667 1.16667 1.66667 (per year length)
Discount to present value for each date : 0.99105 0.96288 0.93401 0.90509
Coupon rate    : 6.00 %
Bond VALUE     : 101.88819
Bond duration  : 1.58080
Bond convexity : 2.59243

Bond yield : 0.05975
Par yield coupon : 0.05004


Process returned 0 (0x0)   execution time : 45.476 s
Press any key to continue.

[Execution time is as large as 45.476 s because of the time needed to manually input starting data into the interactive prompt]

Answer 1

Bond::Bond(valarray<double> inp_times, const double inp_coupon)
    : cash_times{inp_times}, coupon_rate{inp_coupon}, discount_factors(cash_times.size()),
    bond_value{0}, duration{0}, convexity{0}, yield{0}, par_yield_coupon{0}
{}

Why are you passing imp_times by value? If you meant to use the "sink" constructor idiom, you forgot the std::move. If not, you forgot that C++ has value semantics?

This constructor body should go inline in the header, inside the class definition, not in the CPP file.

Use inline initializers for all the members that are not set by arguments. E.g. the class would define the member as: double bond_value = 0.0; inside the class. Then, don't list them all in the constructor init list.

BTW, don't use double for money!!!

---

inputFromUser, with a fist full of "out" parameters, is not very nice. It would be awkward to call too, having to first declare all those things and then pass them to the function, rather than initializing (and maybe making const).

---

Use switch instead of a cascade of if/else when you want to do something based on the value of an enumeration variable.

---

~KnownRateBond() = default;
There's no need to do that in the derived classes. You only needed =default in the base class in order to mark it as virtual. Declaring the destructor -- even as =default -- disables the auto-generated copy constructor and assignment operators.

---

(BTW, I see a lot of good things too!)