• C++ Programming for Financial Engineering
    Highly recommended by thousands of MFE students. Covers essential C++ topics with applications to financial engineering. Learn more Join!
    Python for Finance with Intro to Data Science
    Gain practical understanding of Python to read, understand, and write professional Python code for your first day on the job. Learn more Join!
    An Intuition-Based Options Primer for FE
    Ideal for entry level positions interviews and graduate studies, specializing in options trading arbitrage and options valuation models. Learn more Join!

Binomial Option Pricing Application

Joined
8/5/14
Messages
312
Points
303
Yesterday I cracked upon Hull Chapter 13 and read the chapter on binomial trees. I decided it would be a cool project to code it up. About 6 hours later, I've got a nice Binomial Option (call) pricing application (I could easily extend it for put options later).

What I learned: In the end, the application showed that the Binomial Tree model converged to the exact solution as the number of time steps got large.

How I did it: First I wrote a general binomial lattice template class. Then I wrote a class which contained all the Binomial Tree data (using composition of the lattice) needed in order to solve the problem (i.e. find the option price today). The lattice data member stored a tuple<double, double>, which contained the stock price and option price of the respective position in the tree. Algorithms were implemented to solve the problem.

I was wondering if any of the experts (@APalley , @Daniel Duffy ) could give me some feedback on my implementation in C++. I know I get the right answer, but I would like to improve efficiency. Thanks in advance.

Binomial Lattice Class:
Code:
// Ryan Shrott: Aug 11, 12
#ifndef BinomialLattice_h
#define BinomialLattice_h

#include <vector>
#include <iostream>

template <typename T>
class BinomialLattice
{
private:
    std::vector< std::vector<T> > m_lat;
    unsigned m_timeSteps;
    unsigned m_intial;
public:
    BinomialLattice();
    BinomialLattice(unsigned timeSteps, const T& intial);
    BinomialLattice(const BinomialLattice<T>& rhs);
    virtual ~BinomialLattice();

    // Assignment operator
    BinomialLattice<T>& operator= (const BinomialLattice<T>& rhs);

    // Data fetching (and Modifying, since we return by reference) (Remark: A data pounsigned is uniquely defined by a timeStep and a state)
    T& operator()(const unsigned& timeStep, const unsigned& state); // define the 0th state to be the bottom most state of lattice (for fixed time step)
    const T& operator()(const unsigned& timeStep, const unsigned& state) const; // const method of fetching state data for some time step

    // Accessing size of lattice
    unsigned get_timeSteps() const;

    // Access and change member variables
    const std::vector< std::vector<T> >& Lat(void) const;
    std::vector< std::vector<T> >& Lat(void);

    unsigned Steps(void) const;

    // Mutate member variables
    void Lat(const std::vector< std::vector<T> > lat);
    void Steps(const unsigned steps);

    // Print Lattice data to console window
    template <typename T>
    friend std::ostream& operator<< (std::ostream& os, const BinomialLattice<T>& lat);
};

#ifndef BinomialLattice_cpp
#include "BinomialLattice.cpp"
#endif

#endif

Code:
// Ryan Shrott: Aug 11, 12
#ifndef BinomialLattice_cpp
#define BinomialLattice_cpp

#include "BinomialLattice.h"

template<typename T>
BinomialLattice<T>::BinomialLattice()
{
    m_lat.resize(2 + 1); // Define time(0) = first time step, so size is timeSteps + 1
    for (unsigned i = 0; i < m_lat.size(); ++i)
    {
        m_lat[i].resize(i + 1, 0.0); // The i'th time step has i+1 possible states
    }
    m_timeSteps = 2;
}

template <typename T>
BinomialLattice<T>::BinomialLattice(unsigned timeSteps,const T& intial)
{
    m_lat.resize(timeSteps + 1); // Define time(0) = first time step, so size is timeSteps + 1
    for (unsigned i = 0; i < m_lat.size(); ++i)
    {
        m_lat[i].resize(i + 1, intial); // The i'th time step has i+1 possible states
    }
    m_timeSteps = timeSteps;
}

template <typename T>
BinomialLattice<T>::BinomialLattice(const BinomialLattice<T> & rhs)
{
    m_lat = rhs.m_lat; // vector assignment overload
    m_timeSteps = rhs.m_timeSteps;
}

template <typename T>
BinomialLattice<T>::~BinomialLattice()
{
}

template<typename T>
BinomialLattice<T>& BinomialLattice<T>::operator=(const BinomialLattice<T>& rhs)
{
    if (this != &rhs)
    {
        m_timeSteps = rhs.m_timeSteps;
        m_lat.resize(m_timeSteps + 1); //resize time steps
        for (unsigned i = 0; i < m_lat.size(); ++i)
        {
            m_lat[i].resize(i + 1, 0.0); // The i'th time step has i+1 possible states
        }
        for (unsigned t = 0; t <= m_timeSteps; ++t) // traversing time
        {
            for (unsigned s = 0; s <= t; ++s) // traversing possible states
            {
                (*this)(t, s) = rhs(t, s);
            }
        }
    }
    return *this;
}

template <typename T>
T& BinomialLattice<T>::operator()(const unsigned & timeStep, const unsigned & state)
{
    return m_lat[timeStep][state];
}

template <typename T>
const T& BinomialLattice<T>::operator()(const unsigned & timeStep, const unsigned & state) const
{
    return m_lat[timeStep][state];
}

template <typename T>
unsigned BinomialLattice<T>::get_timeSteps() const
{
    return m_timeSteps;
}

template<typename T>
const std::vector< std::vector<T> >& BinomialLattice<T>::Lat() const
{
    return m_lat;
}

template<typename T>
std::vector<std::vector<T>>& BinomialLattice<T>::Lat(void)
{
    return m_lat;
}

template<typename T>
unsigned BinomialLattice<T>::Steps() const
{
    return m_timeSteps;
}

template<typename T>
void BinomialLattice<T>::Lat(const std::vector<std::vector<T>> lat)
{
    m_lat = lat;
}

template<typename T>
void BinomialLattice<T>::Steps(const unsigned steps)
{
    m_timeSteps = steps;
}

template <typename T>
std::ostream& operator<< (std::ostream& os, const BinomialLattice<T>& lat)
{
    std::cout << "Lattice: " << std::endl;
    for (unsigned t = 0; t <= lat.get_timeSteps(); ++t) // traversing time
    {
        for (unsigned s = 0; s <= t; ++s) // traversing possible states
        {
            std::cout << lat(t, s) << "\t";
        }
        std::cout << std::endl;
    }
    return os;
}

#endif

Binomial Tree Problem Solver Class:
Code:
// Ryan Shrott: Aug 11, 12
#ifndef BinomialTree_h
#define BinomialTree_h

#include "BinomialLattice.h"
#include <algorithm>
#include <cmath>
#include <boost/tuple/tuple.hpp>    // Tuple class
#include <boost/tuple/tuple_io.hpp> // I/O for tuples

typedef boost::tuple<double, double> data; // format: (stock price, option price)

class BinomialTree // contains all the data which conmprise a binomial tree model for the underlying
{
private:
    BinomialLattice<data> m_lat; // each point of lattice stores the stock price and option price, as a tuple
    unsigned m_timeSteps;
    double m_u; // UP multiplier
    double m_d; // DOWN multiplier
    double m_riskFree; // risk free interest rate
    double m_K; // strike price
    double m_delta_T; // length of a time-step
    double m_S_zero; // intitial value of stock
public:
    BinomialTree(double timeExpiration, double numSteps, double sigma, double r, double k, double s);
    ~BinomialTree();
    BinomialTree(const BinomialTree& rhs);
    BinomialTree& operator=(const BinomialTree& rhs);

    // Solve option prices director
    void Solve(void); // calculates all option prices at each time step and modify lattice structure

    // Calculate Risk Free probability measure
    double RiskFreeProb(void) const;

    // Data fetching (and Modifying, since we return by reference) (Remark: A data pounsigned is uniquely defined by a timeStep and a state)
    data& operator()(const unsigned& timeStep, const unsigned& state); // define the 0th state to be the bottom most state of lattice (for fixed time step)
    const data& operator()(const unsigned& timeStep, const unsigned& state) const; // const method of fetching state data for some time step

    // Accessing size of lattice
    unsigned get_timeSteps() const;

    // Print Lattice data to console window
    friend std::ostream& operator<< (std::ostream& os, const BinomialTree& lat);
};

#endif

Code:
// Ryan Shrott: Aug 11, 12
#include "BinomialTree.h"

BinomialTree::BinomialTree(double timeExpiration, double numSteps, double sigma, double r, double k, double s)
{
    m_timeSteps = numSteps;
    m_delta_T = timeExpiration/ (numSteps + 1);
    m_lat = BinomialLattice<data>(m_timeSteps, boost::make_tuple(0.0,0.0));
    m_u = exp(sigma * sqrt(m_delta_T));
    m_d = exp(-1 * sigma * sqrt(m_delta_T));
    m_riskFree = r;
    m_S_zero = s;
    m_K = k;

    // Assigning Stock Prices at each time step
    for (unsigned t = 0; t <= this->get_timeSteps(); ++t) // traversing time
    {
        for (unsigned s = 0; s <= t; ++s) // unsign stock prices as a function of stock price, up multiplier and down multiplier
        {
            (*this)(t, s).get<0>() = m_S_zero * pow(m_u, s) * pow(m_d, t - s);
        }
    }
    // Assigning option prices for last time step
    for (unsigned s = 0; s <= this->get_timeSteps(); ++s)
    {
        (*this)(this->get_timeSteps(), s).get<1>() =( (*this)(this->get_timeSteps(), s).get<0>() > m_K ? (*this)(this->get_timeSteps(), s).get<0>() - m_K: 0.0 ); // intialize option prices at each time step to zero
    }
}

BinomialTree::~BinomialTree()
{
}

BinomialTree::BinomialTree(const BinomialTree & rhs)
{
    m_lat.BinomialLattice<data>::operator=(rhs.m_lat);
    m_u = rhs.m_u;
    m_d = rhs.m_d;
    m_riskFree = rhs.m_riskFree;
    m_delta_T = rhs.m_delta_T;
    m_S_zero = rhs.m_S_zero;
}

BinomialTree & BinomialTree::operator=(const BinomialTree & rhs)
{
    if (this != &rhs)
    {
        m_lat.BinomialLattice<data>::operator=(rhs.m_lat);
        m_u = rhs.m_u;
        m_d = rhs.m_d;
        m_riskFree = rhs.m_riskFree;
        m_delta_T = rhs.m_delta_T;
        m_S_zero = rhs.m_S_zero;
    }
    return *this;
}

void BinomialTree::Solve(void)
{
    for (unsigned t = this->get_timeSteps() - 1; t != -1; --t) // traversing time backwards
    {
        for (unsigned s = 0; s <= t; ++s) // unsign stock prices as a function of stock price, up multiplier and down multiplier
        {
            if ((*this)(t + 1, s + 1).get<1>() == 0 && (*this)(t + 1, s).get<1>() == 0)
            {
                (*this)(t, s).get<1>() = 0.0;
            }
            else // the expected value (using risk free prob) discounted continuously
            {
                (*this)(t, s).get<1>() = exp(-m_riskFree * m_delta_T) *
                                        (this->RiskFreeProb() * (*this)(t + 1, s + 1).get<1>() + (1 - this->RiskFreeProb()) * (*this)(t + 1, s).get<1>());
            }
        }
    }

}

double BinomialTree::RiskFreeProb(void) const
{
    return (exp(m_riskFree * m_delta_T) - m_d) / (m_u - m_d);
}

data& BinomialTree::operator()(const unsigned& timeStep, const unsigned& state)
{
    return m_lat.Lat()[timeStep][state]; // return a double tuple
}

const data& BinomialTree::operator()(const unsigned & timeStep, const unsigned & state) const
{
    return m_lat.Lat()[timeStep][state]; // return a double tuple
}

unsigned BinomialTree::get_timeSteps() const
{
    return m_timeSteps;
}

std::ostream& operator<< (std::ostream& os, const BinomialTree& lat)
{
    std::cout << "Lattice: " << std::endl;
    for (unsigned t = 0; t <= lat.get_timeSteps(); ++t) // traversing time
    {
        for (unsigned s = 0; s <= t; ++s) // traversing possible states
        {
            std::cout << lat(t, s) << "\t";
        }
        std::cout << std::endl;
    }
    return os;
}

Main Function for Testing:
Code:
// Ryan Shrott: Aug 11, 12

#include "BinomialTree.h"

#include <iostream>
using namespace std;
#include <cmath>
#include <boost/tuple/tuple.hpp>    // Tuple class
#include <boost/tuple/tuple_io.hpp> // I/O for tuples

typedef boost::tuple<double, double> data; // data = (stock price, option price)

int main()
{
    double numSteps, timeToExp;
    cout << "Enter time to expiration (T): "; cin >> timeToExp;
    cout << "Enter the number of time steps: "; cin >> numSteps;
    BinomialTree TreeProblem(timeToExp, numSteps, 0.30, 0.08, 65.0, 60.0); // notation: BinomialTreedouble(timeExpiration, numSteps, sigma, r, K, S_0)
    TreeProblem.Solve();
    //cout << TreeProblem << endl;
    cout << "The option price today should be: " << TreeProblem(0, 0).get<1>() << endl;
    return 0;
}
 
A few comments after a brief glance:
- `#include <cmath>` (as opposed to `#include <math.h>`) followed `exp` is an error: use `std::exp` instead // same for `sqrt`, `pow`
- `#include "BinomialLattice.cpp"` -- you don't normally want to do that in a header file, unless this code was meant to use the Boost's header-only `ipp` template interface-implementation separation pattern (but then it should have explicitly used the `ipp` in the name)
- I'd have used `std::size_t` instead of `unsigned` -- as the standard size type it expresses the intent better
- you never want to use `const unsigned&` for an argument (use either `unsigned` or `const unsigned` instead, there's more to say on the latter variant); if this was meant to be an "optimization", it's worth understanding why this doesn't make sense -- pointer itself (references compile to const-pointers) will be 64-bit on the current platforms (whether you're targeting LLP64/Windows or LP64/POSIX), while `unsigned` for (L)LP64 is 32-bit -- more importantly this adds an extra indirection and potentially even prevents compiler optimization due to the introduced aliasing; in general, I'd avoid using language features (even references) for "optimization" before understanding the reasons _why_ -- in this case thinking through the _why_ step would prevent resorting to it in the first place
- `typedef boost::tuple<double, double> data` combined with the necessity of documenting with a comnent `data = (stock price, option price)` indicates that you should probably use a `struct` here (with the corresponding two members): if you have to document the usage intent in a comment, this implies that the chosen data type is too poor to express the intent itself; alternatively, making `stock_price` and `option_price` types themselves and using C++14 type-addressed tuple could be another thing to explore -- C++14 - Wikipedia, the free encyclopedia -- my initial instinct would be to start with a `struct`, though
- related to the aforementioned choice of using a `tuple`: `(*this)(this->get_timeSteps(), s).get<1>() =( (*this)(this->get_timeSteps(), s).get<0>() > m_K ? (*this)(this->get_timeSteps(), s).get<0>() - m_K: 0.0 ); // intialize option prices at each time step to zero` -- are you fully comfortable with the clarity of this line of code?
- `get_timeSteps` doesn't fit any of the mainstream naming/capitalization conventions (like `getTimeSteps` or `get_time_steps` would)
- `void Solve(void)` -- this is C-ism; in C++ we write `void Solve()`
- `virtual ~BinomialLattice();` -- this is a publicly announced commitment that your API design requires that `BinomialLattice` will be used as a parent class (i.e., the intent that it will be inherited-from); where are the derived classes? I can see that BinomialTree HAS-A BinomialLattice (as opposed to saying that `BinomialTree` IS-A `BinomialLattice`; and, good, composition is usually preferred to inheritance), so this seems to be a contradiction in the design.
- BinomialTree TreeProblem(timeToExp, numSteps, 0.30, 0.08, 65.0, 60.0); // notation: BinomialTreedouble(timeExpiration, numSteps, sigma, r, K, S_0)` -- is this interface "Easy to Use Correctly and Hard to Use Incorrectly"?
 
Last edited:
- `#include "BinomialLattice.cpp"` -- you don't normally want to do that in a header file, unless this code was meant to use the Boost's header-only `ipp` template interface-implementation separation pattern (but then it should have explicitly used the `ipp` in the name)

I do this because BinomialLattice class is a template. So instead of including the .cpp file in main, I just include it in the .h file and then I can include the .h file in main().
 
- you never want to use `const unsigned&` for an argument (use either `unsigned` or `const unsigned` instead, there's more to say on the latter variant); if this was meant to be an "optimization", it's worth understanding why this doesn't make sense -- pointer itself (references compile to const-pointers) will be 64-bit on the current platforms (whether you're targeting LLP64/Windows or LP64/POSIX), while `unsigned` for (L)LP64 is 32-bit -- more importantly this adds an extra indirection and potentially even prevents compiler optimization due to the introduced aliasing; in general, I'd avoid using language features (even references) for "optimization" before understanding the reasons _why_ -- in this case thinking through the _why_ step would prevent resorting to it in the first place

Interesting. my thought of doing this was that const references were always more efficient that call by value.
 
- `virtual ~BinomialLattice();` -- this is a publicly announced commitment that your API design requires that `BinomialLattice` will be used as a parent class (i.e., the intent that it will be inherited-from); where are the derived classes? I can see that BinomialTree HAS-A BinomialLattice (as opposed to saying that `BinomialTree` IS-A `BinomialLattice`; and, good, composition is usually preferred to inheritance), so this seems to be a contradiction in the design.

Shouldn't be virtual here. And yes, I used composition here, not inheritance.
 
I do this because BinomialLattice class is a template. So instead of including the .cpp file in main, I just include it in the .h file and then I can include the .h file in main().

Yeah, that's the Boost's `.ipp` files:
In the C++ Boost libraries, why is there a ".ipp" extension on some header files
Boost - Users - ipp files

I'd also generally reserve `.h` for C (not C++ -- use `.hpp` for that); although, sadly, not everyone uses this convention, there are good reasons to do so (including that C and C++ are fundamentally different languages, we wouldn't use `.c` instead of `.cpp`, either):
*.h or *.hpp for your class definitions

Interesting. my thought of doing this was that const references were always more efficient that call by value.
Yeah, hence the general advice: "I'd avoid using language features (even references) for "optimization" before understanding the reasons _why_." In this context: Try to think about the reasons why passing-by-reference would be more efficient than passing-by-value -- this should lead to the answer to when it's more efficient (hint: not always)... and when it's not :)
// Note the part about "big" objects: Standard C++
There's one other piece of advice from Andrei Alexandrescu:
  • The only good intuition: “I should time this.”
Three Optimization Tips for C++—Andrei Alexandrescu : Standard C++

Overall, this is a good exercise, starting with a piece of code and improving; if you're interested, you may also get more comments here:
Newest 'c++' Questions
Although generally C++11 and newer will be preferred -- this, however, gives you an opportunity to also try to practice that -- e.g., `using ` instead of `typedef` (it's much more readable anyway): Type alias, alias template (since C++11) - cppreference.com
 
"I've got a nice Binomial Option (call) pricing application (I could easily extend it for put options later)."

Some remarks/idea on the problem:

1. How do you extend to put options? e.g. using an OO or functional approach?
2. The class is hard-coded for GBM and equities.
3. The alias template is useful (better than typedef).
4. What about a lattice structure ADT that subsumes both binomial and trinomial? If you are ambitious try variadic tuples :)
5. I would like to use the code ADT as data to do polynomial (e.g. Neville) interpolation.
6. Using lambda functions to configure the application (e.g. payoffs).
7. How to support early exercise,dividends and barriers.
 
Last edited:
I do this because BinomialLattice class is a template. So instead of including the .cpp file in main, I just include it in the .h file and then I can include the .h file in main().
One option is that .hpp includes .cpp at the former's end and then main() just use .hpp.

I don't like .ipp. Anyways, each to his own.
 
"The only good intuition: “I should time this.”"

Agree 90%. The word 'only' sticks in my craw. It sounds like ad hoc.
One can use science to predict the efficiency.
 
Last edited:
I'd also generally reserve `.h` for C (not C++ -- use `.hpp` for that); although, sadly, not everyone uses this convention, there are good reasons to do so (including that C and C++ are fundamentally different languages, we wouldn't use `.c` instead of `.cpp`, either):

The only reason I do that is because Visual Studio 2015 defaults to .h when I create a new file: I just type in the file name, and the outcome is .h.
 
1. How do you extend to put options? e.g. using an OO or functional approach?

I think it's easiest to just have a SolveCall(), and SolvPut() function, which is what I did. There's a lot of repeating code, but it gets the job done.

4. What about a lattice structure ADT that subsumes both binomial and trinomial? If you are ambitious try variadic tuples :)

But my code already produces very accurate results for many batch trials. So why is that needed?
5. I would like to use the code ADT as data to do polynomial (e.g. Neville) interpolation.

What's to interpolate?

7. How to support early exercise,dividends and barriers.

I'll think about this :)
 
"But my code already produces very accurate results for many batch trials. So why is that needed?"
Famous last words. :) quants spend 2 days writing code and 3 months debugging it.

You tested it a few benign test cases?

Try e.g. small vol and large drift etc.
//
If your lattice is generic you can do all these with no change of code
Neville's algorithm - Wikipedia, the free encyclopedia
Romberg's method - Wikipedia, the free encyclopedia
Pascal's triangle - Wikipedia, the free encyclopedia

Cool?
 
Last edited:
" There's a lot of repeating code, but it gets the job done. "
Copy and past syndrome? Yikes :)

See how I resolved put and call in my code (they are just functions, so you don't need call(), put().

BackwardInduction() takes an arbitrary function at T, and you configure in main() with a lambda etc. chap 11, page 8, page 9.
 
"I know I get the right answer, but I would like to improve efficiency."

Very good but only (less than) half the battle. 90% is maintaining software as it is extended, debugged, falls apart down all the years etc.

Some trader will ask you to add stuff to it.

Or like it works for most vol but not when vol = 62.13. Simple?
 
Last edited:
But my code already produces very accurate results for many batch trials. So why is that needed?

One more idea: How accurate? What are the expected error bounds given the steps count (and the market parameters, model parameters, contract parameters)?

You can document this via an executable specification which can be run continuously when you make changes (new products/contracts, new algorithms, different markets) -- i.e., a test suite.
I'd recommend Catch: philsquared/Catch · GitHub
Here's how to get started: Catch/tutorial.md at master · philsquared/Catch · GitHub
I'd recommend using BDD style: Catch/tutorial.md at master · philsquared/Catch · GitHub
Testdriven C++ with Catch - Phil Nash:

// Possible next step: property-based testing with RapidCheck (integrates with Catch): Generating test cases so you don’t have to

If you have never done testing before, you can think of this as "comments on steroids" -- in that tests are executable and stay in sync with the code (and can be randomly generated to test properties beyond individual cases, see RapidCheck).
IME the "hierarchy of needs" can be described as follows:
- Ideally, if you can express something in a type (like `std::size_t` instead of `int`), that's the best way to assure correctness at compile-time. It goes far beyond that and is one of the major advantages of statically typed languages; see "Type Driven Design" and "Why type-first development matters".
- Then, tests to assure correctness at development / pre-shipping run-time.
- Then, asserts / internal checks / exceptions for deployment / post-shipping run-time.
- Comments, docs -- probably won't be read, probably won't be maintained, most likely won't stay in sync with the code // well, OK, it gets better sometimes ;-)
 
Last edited:
"One more idea: How accurate? What are the expected error bounds given the steps count (and the market parameters, model parameters, contract parameters)?"

Depends on if CRR, JR, Tian, Leisin-Reimer etc. See Stefanie Mueller's thesis.
(CRR is very bad). Do you want greeks?

Binomial is 1st order accurate at best; trinomial is better.

http://webdoc.sub.gwdg.de/ebook/dissts/Kaiserslautern/Mueller2009.pdf
 
Last edited:
Back
Top