// CRR_Model.cpp
//
// Cox-Ross-Rubinstein Binomial Option Pricing Model
//
// Supports:
// - European and American options
// - Calls and puts
// - Stock and Futures underlyings
//
// Mathematical background:
//
// Stock tree:
//
// S(i,j) = S0 * u^j * d^(i-j)
//
// where:
//
// u = exp(sigma * sqrt(dt))
// d = 1/u
//
// Risk neutral probability:
//
// p = ( exp((r-q)dt) - d ) / (u-d)
//
// Discount factor:
//
// disc = exp(-r*dt)
//
// Futures options:
//
// Futures prices are martingales under the risk-neutral measure.
// Therefore the futures tree is obtained by:
//
// F(i,j) = p * F(i+1,j+1) + (1-p) * F(i+1,j)
//
// with no discounting.
//
// The option itself is still discounted using the risk-free rate.
//
#include <iostream>
#include <vector>
#include <cmath>
#include <algorithm>
#include <iomanip>
#include <string>
// ------------------------------------------------------------
// Enumerations
// ------------------------------------------------------------
enum Exercise { EUROPEAN, AMERICAN };
enum OptionType { CALL, PUT };
enum Underlying { STOCK, FUTURE };
// ------------------------------------------------------------
// CRR Model Class
// ------------------------------------------------------------
class CRR_Model
{
public:
CRR_Model( double S, double K, double sigma, double r, double q, double T, Exercise exercise,
OptionType optionType, Underlying underlying, int n, int m = 0)
:
S_(S), K_(K), sigma_(sigma), r_(r), q_(q), T_(T),
exercise_(exercise), optionType_(optionType), underlying_(underlying),
n_(n), m_(m)
{
// Time step is determined by the underlying tree.
//
// For stocks: option maturity = underlying maturity
// For futures: futures expiry determines the lattice spacing
int periods = (underlying_ == Underlying::FUTURE) ? m_ : n_;
dt_ = T_ / periods;
u_ = std::exp( sigma_ * std::sqrt(dt_) );
d_ = 1/u_;
disc_ = std::exp(-r_*dt_);
// Risk neutral probability
p_ = ( std::exp( (r_-q_)*dt_)-d_ ) / (u_-d_);
earliestExercisePeriod_=-1;
}
// --------------------------------------------------------
// Price the option
// --------------------------------------------------------
double price()
{
earliestExercisePeriod_=-1;
int depth = (underlying_==Underlying::STOCK) ? n_ : m_;
auto stockTree = buildStockTree(depth);
std::vector<std::vector<double>> underlyingTree;
if(underlying_ == Underlying::STOCK)
{
underlyingTree = stockTree;
}
else
{
underlyingTree = buildFuturesTree(stockTree);
}
// Terminal payoff
std::vector<double> V(n_+1);
for(int j = 0; j <= n_; ++j)
{
V[j]=intrinsic( underlyingTree[n_][j] );
}
// Backward induction
const double exerciseTolerance = 1e-12;
for(int i = n_-1; i >= 0; --i)
{
for(int j = 0; j <= i; ++j)
{
double continuation = disc_ * ( p_*V[j+1] + (1-p_)*V[j] );
if(exercise_ == Exercise::AMERICAN)
{
double exerciseValue = intrinsic( underlyingTree[i][j] );
if( exerciseValue > continuation + exerciseTolerance )
{
earliestExercisePeriod_=i;
}
V[j] = std::max( continuation, exerciseValue );
}
else
{
V[j] = continuation;
}
}
}
return V[0];
}
int optimalExercise()
{
return earliestExercisePeriod_;
}
void printResults(std::string name)
{
std::cout << std::left << std::setw(45) << name << price() << "\n";
}
private:
// Input parameters
double S_; // Initial stock/futures price
double K_; // Strike price
double sigma_; // Volatility
double r_; // Risk-free rate
double q_; // Dividend yield
double T_; // Time to expiry
Exercise exercise_;
OptionType optionType_;
Underlying underlying_;
int n_; // Option periods
int m_; // Futures periods
// Model quantities
double dt_;
double u_;
double d_;
double disc_;
double p_;
int earliestExercisePeriod_;
// Build stock price tree
std::vector<std::vector<double>> buildStockTree(int periods)
{
std::vector<std::vector<double>> tree(periods + 1);
for(int i = 0; i <= periods; ++i)
{
tree[i].resize(i+1);
for(int j = 0; j <= i; ++j)
{
tree[i][j] = S_ * std::pow(u_,j) * std::pow(d_, i-j);
}
}
return tree;
}
// Construct futures lattice
//
// At futures expiry:
//
// F(T,T)=S(T)
//
// We then roll backwards because futures prices
// are martingales.
std::vector<std::vector<double> >
buildFuturesTree(std::vector<std::vector<double> >& stockTree)
{
std::vector<std::vector<double> > futureTree(m_+1);
for(int i = 0; i <= m_; ++i)
futureTree[i].resize(i+1);
// Terminal futures prices equal stock prices
for(int j = 0; j <= m_; ++j)
futureTree[m_][j]=stockTree[m_][j];
// Futures martingale rollback
for(int i = m_-1; i >= 0; --i)
{
for(int j = 0; j <= i; ++j)
{
futureTree[i][j] = p_ * futureTree[i+1][j+1] + (1 - p_) * futureTree[i+1][j];
}
}
return futureTree;
}
double intrinsic(double price)
{
if(optionType_ == OptionType::CALL)
return std::max( price - K_, 0.0 );
return std::max( K_- price, 0.0 );
}
};
// ------------------------------------------------------------
// Test Cases
// ------------------------------------------------------------
int main()
{
double S=100;
double K=110;
double sigma=0.30;
double r=0.02;
double q=0.01;
double T=0.25;
std::cout << std::fixed << std::setprecision(4);
std::cout << "--- Equity Options ---\n";
CRR_Model europeanCall( S, K, sigma, r, q, T, EUROPEAN, CALL, STOCK, 15 );
europeanCall.printResults("European Call");
CRR_Model europeanPut( S, K, sigma, r, q, T, EUROPEAN, PUT, STOCK, 15 );
europeanPut.printResults("European Put");
CRR_Model americanCall( S, K, sigma, r, q, T, AMERICAN, CALL, STOCK, 15 );
americanCall.printResults( "American Call (expected ~2.60)" );
CRR_Model americanPut( S, K, sigma, r, q, T, AMERICAN, PUT, STOCK, 15 );
americanPut.printResults( "American Put (expected ~12.36)" );
std::cout << "\n--- Futures Options ---\n";
CRR_Model futuresCall( S, K, sigma, r, q, T, EUROPEAN, CALL, FUTURE, 10, 15 );
futuresCall.printResults( "European Futures Call" );
CRR_Model americanFutureCall( S, K, sigma, r, q, T, AMERICAN, CALL, FUTURE, 10, 15 );
americanFutureCall.printResults( "American Futures Call (expected ~1.66)" );
std::cout << "Optimal exercise period = " << americanFutureCall.optimalExercise() << " (expected 7)\n";
return 0;
}
