#include <iostream>
#include <cmath>
#include <utility>
#include <vector>
#include <cassert>
#include <memory>

using namespace std;

// A Portfolio indicates a certain amount of some underlying product (e.g. Microsoft stock), and a certain amount of
// cash (e.g. US dollars).
struct Portfolio {
    // The amount of underlying we are long or short (this is the traditional "delta" of our option)
    double underlying;
    // The amount of currency we should be long or short
    double currency;

    // Calculate the currencyValue if 1 underlying = spot * currency
    double currencyValue(double spot) const {
        return underlying * spot + currency;
    }
};

ostream& operator<<(ostream& os, const Portfolio &portfolio) {
    return os << "(underlying=" << portfolio.underlying << " currency=" << portfolio.currency << ")";
}

//
// A Rates object indicates how a portfolio grows over time.
//
class Rates {
    double period_;
    double underlyingRate_;
    double currencyRate_;

    double discount(double amount, double rate) const {
        return exp(-rate * period_) * amount;
    }

    public:
    Rates(double period, double underlyingRate, double currencyRate)
        : period_(period),
        underlyingRate_(underlyingRate),
        currencyRate_(currencyRate) {
    }

    Portfolio discount(const Portfolio& portfolio) const {
        Portfolio result;
        result.underlying = discount(portfolio.underlying, underlyingRate_);
        result.currency = discount(portfolio.currency, currencyRate_);
        return result;
    }
};

// A spot process indicates how much a spot price can move up or down.
struct SpotProcess {
    virtual double calcSpotUp(double spot) const = 0;
    virtual double calcSpotDown(double spot) const = 0;
};

// A GeometricSpotProcess allows the spot price to move up or down by a certain scale factor.
class GeometricSpotProcess : public SpotProcess {
    double scaleFactor_;

    public:
    GeometricSpotProcess(double period, double sigma) 
        : scaleFactor_(exp(sqrt(period) * sigma)) {
    }

    double calcSpotUp(double spot) const {
        return spot * scaleFactor_;
    }

    double calcSpotDown(double spot) const {
        return spot / scaleFactor_;
    }
};

// ExpiryFunction determines what is the required payoff at expiration as a function of spot price at expiration.
struct ExpiryFunction {
    virtual Portfolio getExpiryPortfolio(double spot) const = 0;
};

struct BinomialSolver {
    virtual Portfolio solve(int steps, double spot) const = 0;
    virtual ~BinomialSolver() {}
};

struct CallOption : public ExpiryFunction {
    const double strike_;

    public:
    CallOption(double strike)
        : strike_(strike) {
    }

    // ExpiryFunction
    Portfolio getExpiryPortfolio(double spot) const {
        Portfolio result;

        // If the stock is worth more than what we can pay for it due to our option, we're going to exchange the underlying for the stock.
        // Otherwise, we're not going to do anything, and the option is going to lapse.
        if(spot > strike_) {
            // We obtain 1 underlying
            result.underlying = 1.0;
            // We pay out strike price amount of dollars
            result.currency = -strike_;
        } else {
            // We'd lose money if we exercised the option, so we just let it expire, and have nothing.
            result.underlying = 0.0;
            result.currency = 0.0;
        }
        return result;
    }
};

// The most "direct" approach to solving this problem is to perform it recursively.
class RecursiveSolver : public BinomialSolver {
    const Rates rates_;
    const SpotProcess& spotProcess_;
    const ExpiryFunction& expiryFunction_;

    public:
    RecursiveSolver(const Rates& rates, const SpotProcess& spotProcess, const ExpiryFunction& expiryFunction)
        : rates_(rates),
        spotProcess_(spotProcess),
        expiryFunction_(expiryFunction) {
    }

    // Determine what equivalent replicating portfolio will have exactly the same payoff as the expiry function at expiration
    // Our strategy is to have that equivalent replicating portfolio, and each time the spot price moves, we exchange
    // one product for the other at the new spot price in a certain amount, such that once we are at expiration, no matter
    // which path the spot has moved in, our replicating portfolio has exactly the same payoff as the expiry function.
    //
    // steps: number of steps to expiration
    // spot: current spot price
    // returns: replicating portfolio that we should have to be able to re-hedge at zero-cost right through to expiration, and get the expiry portfolio.
    Portfolio solve(int steps, double spot) const {
        if(steps == 0)
            return expiryFunction_.getExpiryPortfolio(spot);


        double spotUp = spotProcess_.calcSpotUp(spot);
        double spotDown = spotProcess_.calcSpotDown(spot);

        double currencyValueUp, currencyValueDown;

        // What currencyValue must our replicating portfolio have if the spot moves up/down? 

        // This is horrendously inefficient if our spotProcess does in fact re-combine:
        // (i.e. spotProcess.calcSpotDown(spotProcess.calcSpotUp(s)) == spotProcess.calcSpotUp(spotProcess.calcSpotDown(s))).
        // Does this mean that this model correctly handles non-combining spot processes? I dunno...
        // We could optimize to take advantage of combining spot processes by keeping a cache keyed off (steps, spot)
        // and the required replicating portfolio for that step. Such a cache would make our algorithm have performance of
        // O(steps ^ 2) _if_ the spotProcess recombines. But the point of this exercise is genericity, not efficiency.
        currencyValueUp = solve(steps - 1, spotUp).currencyValue(spotUp);
        currencyValueDown = solve(steps - 1, spotDown).currencyValue(spotDown);

        // Work out the portfolio we need to have in the future which will have the same currencyValue
        // as the currencyValueUp if the spot price goes up, and the same currencyValue as currencyValueDown if the spot
        // price goes down.

        // That is, take the following two simultaneous equations:
        // futureHedge.underlying * spotUp + futureHedge.currency = currencyValueUp
        // futureHedge.underlying * spotDown + futureHedge.currency = currencyValueDown
        // and re-arrange to get futureHedge.underlying and futureHedge.currency on the left-hand side, and
        // you'll get the following:
        Portfolio futureHedge;
        futureHedge.underlying = (currencyValueUp - currencyValueDown) / (spotUp - spotDown);
        futureHedge.currency = currencyValueUp - futureHedge.underlying * spotUp;
        
        // Now, use rates to figure out how much that portfolio is worth today:
        return rates_.discount(futureHedge);
    }
};

void main() {
    // Evaluate 20 steps
    int steps = 20;

    // 1 year to expiration
    double timeToExpiry = 1.0;

    // Sigma = volatility = 20%
    double sigma = 0.2;

    // Current spot price - is always 100.0 ;-)
    double spot = 80.0;

    // This solver is so inefficient (it's an exponential time algorithm: O(2 ^ steps)),
    // you'll probably have to just kill the program before it reaches 30 iterations... sorry ;-)
    for(int steps = 0; steps != 30; ++steps) {
        double period = timeToExpiry / steps;
        GeometricSpotProcess spotProcess(period, sigma);
        // Let's say that we can borrow the underlying at a rate of 1%, and the currency at a rate of 5%
        Rates rates(period, 0.01, 0.05);
        CallOption callOption(100.0);

        RecursiveSolver solver(rates, spotProcess, callOption);
        Portfolio result = solver.solve(steps, spot);

        // We now know what the equivalent replicating portfolio is. Assuming we have gotten our spot process right (i.e. we have the
        // correct volatility), the currencyValue of this replicating portfolio is the same as the currencyValue
        // of the option (otherwise, we would buy one, and sell the other, and make money).
        cout << "steps = " << steps << " spot = " << spot << " result = " << result << " result currencyValue = " << result.currencyValue(spot) << endl;
    }
}