### A free/open-source library for quantitative finance

DiscreteHedging.cpp

This is an example of using the QuantLib Monte Carlo framework.

    /* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */

    /*  This example computes profit and loss of a discrete interval hedging
        strategy and compares with the results of Derman & Kamal's (Goldman Sachs
        Equity Derivatives Research) Research Note: "When You Cannot Hedge
        Continuously: The Corrections to Black-Scholes"
        http://www.ederman.com/emanuelderman/GSQSpapers/when_you_cannot_hedge.pdf

        Suppose an option hedger sells an European option and receives the
        Black-Scholes value as the options premium.
        Then he follows a Black-Scholes hedging strategy, rehedging at discrete,
        evenly spaced time intervals as the underlying stock changes. At
        expiration, the hedger delivers the option payoff to the option holder,
        and unwinds the hedge. We are interested in understanding the final
        profit or loss of this strategy.

        If the hedger had followed the exact Black-Scholes replication strategy,
        re-hedging continuously as the underlying stock evolved towards its final
        value at expiration, then, no matter what path the stock took, the final
        P&L would be exactly zero. When the replication strategy deviates from
        the exact Black-Scholes method, the final P&L may deviate from zero. This
        deviation is called the replication error. When the hedger rebalances at
        discrete rather than continuous intervals, the hedge is imperfect and the
        replication is inexact. The more often hedging occurs, the smaller the
        replication error.

        We examine the range of possibilities, computing the replication error.
    */

    // the only header you need to use QuantLib
    #include <ql/quantlib.hpp>

    #ifdef BOOST_MSVC
    /* Uncomment the following lines to unmask floating-point
       exceptions. Warning: unpredictable results can arise...

       See http://www.wilmott.com/messageview.cfm?catid=10&threadid=9481
       Is there anyone with a definitive word about this?
    */
    // #include <float.h>
    // namespace { unsigned int u = _controlfp(_EM_INEXACT, _MCW_EM); }
    #endif

    #include <boost/timer.hpp>
    #include <iostream>
    #include <iomanip>

    using namespace QuantLib;

    #if defined(QL_ENABLE_SESSIONS)
    namespace QuantLib {

        Integer sessionId() { return 0; }

    }
    #endif

    /* The ReplicationError class carries out Monte Carlo simulations to evaluate
       the outcome (the replication error) of the discrete hedging strategy over
       different, randomly generated scenarios of future stock price evolution.
    */
    class ReplicationError
    {
    public:
        ReplicationError(Option::Type type,
                         Time maturity,
                         Real strike,
                         Real s0,
                         Volatility sigma,
                         Rate r)
        : maturity_(maturity), payoff_(type, strike), s0_(s0),
          sigma_(sigma), r_(r) {

            // value of the option
            DiscountFactor rDiscount = std::exp(-r_*maturity_);
            DiscountFactor qDiscount = 1.0;
            Real forward = s0_*qDiscount/rDiscount;
            Real stdDev = std::sqrt(sigma_*sigma_*maturity_);
            boost::shared_ptr<StrikedTypePayoff> payoff(
                                                 new PlainVanillaPayoff(payoff_));
            BlackCalculator black(payoff,forward,stdDev,rDiscount);
            std::cout << "Option value: " << black.value() << std::endl;

            // store option's vega, since Derman and Kamal's formula needs it
            vega_ = black.vega(maturity_);

            std::cout << std::endl;

            std::cout << std::setw(8) << " " << " | "
                      << std::setw(8) << " " << " | "
                      << std::setw(8) << "P&L" << " | "
                      << std::setw(8) << "P&L" << " | "
                      << std::setw(12) << "Derman&Kamal" << " | "
                      << std::setw(8) << "P&L" << " | "
                      << std::setw(8) << "P&L" << std::endl;

            std::cout << std::setw(8) << "samples" << " | "
                      << std::setw(8) << "trades" << " | "
                      << std::setw(8) << "mean" << " | "
                      << std::setw(8) << "std.dev." << " | "
                      << std::setw(12) << "formula" << " | "
                      << std::setw(8) << "skewness" << " | "
                      << std::setw(8) << "kurtosis" << std::endl;

            std::cout << std::string(78, '-') << std::endl;
        }

        // the actual replication error computation
        void compute(Size nTimeSteps, Size nSamples);
    private:
        Time maturity_;
        PlainVanillaPayoff payoff_;
        Real s0_;
        Volatility sigma_;
        Rate r_;
        Real vega_;
    };

    // The key for the MonteCarlo simulation is to have a PathPricer that
    // implements a value(const Path& path) method.
    // This method prices the portfolio for each Path of the random variable
    class ReplicationPathPricer : public PathPricer<Path> {
      public:
        // real constructor
        ReplicationPathPricer(Option::Type type,
                              Real strike,
                              Rate r,
                              Time maturity,
                              Volatility sigma)
        : type_(type), strike_(strike),
          r_(r), maturity_(maturity), sigma_(sigma) {
            QL_REQUIRE(strike_ > 0.0, "strike must be positive");
            QL_REQUIRE(r_ >= 0.0,
                       "risk free rate (r) must be positive or zero");
            QL_REQUIRE(maturity_ > 0.0, "maturity must be positive");
            QL_REQUIRE(sigma_ >= 0.0,
                       "volatility (sigma) must be positive or zero");

        }
        // The value() method encapsulates the pricing code
        Real operator()(const Path& path) const;

      private:
        Option::Type type_;
        Real strike_;
        Rate r_;
        Time maturity_;
        Volatility sigma_;
    };

    // Compute Replication Error as in the Derman and Kamal's research note
    int main(int, char* []) {

        try {

            boost::timer timer;
            std::cout << std::endl;

            Time maturity = 1.0/12.0;   // 1 month
            Real strike = 100;
            Real underlying = 100;
            Volatility volatility = 0.20; // 20%
            Rate riskFreeRate = 0.05; // 5%
            ReplicationError rp(Option::Call, maturity, strike, underlying,
                    volatility, riskFreeRate);

            Size scenarios = 50000;
            Size hedgesNum;

            hedgesNum = 21;
            rp.compute(hedgesNum, scenarios);

            hedgesNum = 84;
            rp.compute(hedgesNum, scenarios);

            Real seconds = timer.elapsed();
            Integer hours = int(seconds/3600);
            seconds -= hours * 3600;
            Integer minutes = int(seconds/60);
            seconds -= minutes * 60;
            std::cout << " \nRun completed in ";
            if (hours > 0)
                std::cout << hours << " h ";
            if (hours > 0 || minutes > 0)
                std::cout << minutes << " m ";
            std::cout << std::fixed << std::setprecision(0)
                      << seconds << " s\n" << std::endl;

            return 0;
        } catch (std::exception& e) {
            std::cerr << e.what() << std::endl;
            return 1;
        } catch (...) {
            std::cerr << "unknown error" << std::endl;
            return 1;
        }
    }

    /* The actual computation of the Profit&Loss for each single path.

       In each scenario N rehedging trades spaced evenly in time over
       the life of the option are carried out, using the Black-Scholes
       hedge ratio.
    */
    Real ReplicationPathPricer::operator()(const Path& path) const {

        Size n = path.length()-1;
        QL_REQUIRE(n>0, "the path cannot be empty");

        // discrete hedging interval
        Time dt = maturity_/n;

        // For simplicity, we assume the stock pays no dividends.
        Rate stockDividendYield = 0.0;

        // let's start
        Time t = 0;

        // stock value at t=0
        Real stock = path.front();

        // money account at t=0
        Real money_account = 0.0;

        /************************/
        /*** the initial deal ***/
        /************************/
        // option fair price (Black-Scholes) at t=0
        DiscountFactor rDiscount = std::exp(-r_*maturity_);
        DiscountFactor qDiscount = std::exp(-stockDividendYield*maturity_);
        Real forward = stock*qDiscount/rDiscount;
        Real stdDev = std::sqrt(sigma_*sigma_*maturity_);
        boost::shared_ptr<StrikedTypePayoff> payoff(
                                           new PlainVanillaPayoff(type_,strike_));
        BlackCalculator black(payoff,forward,stdDev,rDiscount);
        // sell the option, cash in its premium
        money_account += black.value();
        // compute delta
        Real delta = black.delta(stock);
        // delta-hedge the option buying stock
        Real stockAmount = delta;
        money_account -= stockAmount*stock;

        /**********************************/
        /*** hedging during option life ***/
        /**********************************/
        for (Size step = 0; step < n-1; step++){

            // time flows
            t += dt;

            // accruing on the money account
            money_account *= std::exp( r_*dt );

            // stock growth:
            stock = path[step+1];

            // recalculate option value at the current stock value,
            // and the current time to maturity
            rDiscount = std::exp(-r_*(maturity_-t));
            qDiscount = std::exp(-stockDividendYield*(maturity_-t));
            forward = stock*qDiscount/rDiscount;
            stdDev = std::sqrt(sigma_*sigma_*(maturity_-t));
            BlackCalculator black(payoff,forward,stdDev,rDiscount);

            // recalculate delta
            delta = black.delta(stock);

            // re-hedging
            money_account -= (delta - stockAmount)*stock;
            stockAmount = delta;
        }

        /*************************/
        /*** option expiration ***/
        /*************************/
        // last accrual on my money account
        money_account *= std::exp( r_*dt );
        // last stock growth
        stock = path[n];

        // the hedger delivers the option payoff to the option holder
        Real optionPayoff = PlainVanillaPayoff(type_, strike_)(stock);
        money_account -= optionPayoff;

        // and unwinds the hedge selling his stock position
        money_account += stockAmount*stock;

        // final Profit&Loss
        return money_account;
    }

    // The computation over nSamples paths of the P&L distribution
    void ReplicationError::compute(Size nTimeSteps, Size nSamples)
    {
        QL_REQUIRE(nTimeSteps>0, "the number of steps must be > 0");

        // hedging interval
        // Time tau = maturity_ / nTimeSteps;

        /* Black-Scholes framework: the underlying stock price evolves
           lognormally with a fixed known volatility that stays constant
           throughout time.
        */
        Calendar calendar = TARGET();
        Date today = Date::todaysDate();
        DayCounter dayCount = Actual365Fixed();
        Handle<Quote> stateVariable(
                              boost::shared_ptr<Quote>(new SimpleQuote(s0_)));
        Handle<YieldTermStructure> riskFreeRate(
                              boost::shared_ptr<YieldTermStructure>(
                                          new FlatForward(today, r_, dayCount)));
        Handle<YieldTermStructure> dividendYield(
                              boost::shared_ptr<YieldTermStructure>(
                                          new FlatForward(today, 0.0, dayCount)));
        Handle<BlackVolTermStructure> volatility(
                              boost::shared_ptr<BlackVolTermStructure>(
                                   new BlackConstantVol(today, calendar, sigma_, dayCount)));
        boost::shared_ptr<StochasticProcess1D> diffusion(
                       new BlackScholesMertonProcess(stateVariable, dividendYield,
                                                     riskFreeRate, volatility));

        // Black Scholes equation rules the path generator:
        // at each step the log of the stock
        // will have drift and sigma^2 variance
        PseudoRandom::rsg_type rsg =
            PseudoRandom::make_sequence_generator(nTimeSteps, 0);

        bool brownianBridge = false;

        typedef SingleVariate<PseudoRandom>::path_generator_type generator_type;
        boost::shared_ptr<generator_type> myPathGenerator(new
            generator_type(diffusion, maturity_, nTimeSteps,
                           rsg, brownianBridge));

        // The replication strategy's Profit&Loss is computed for each path
        // of the stock. The path pricer knows how to price a path using its
        // value() method
        boost::shared_ptr<PathPricer<Path> > myPathPricer(new
            ReplicationPathPricer(payoff_.optionType(), payoff_.strike(),
                                  r_, maturity_, sigma_));

        // a statistics accumulator for the path-dependant Profit&Loss values
        Statistics statisticsAccumulator;

        // The Monte Carlo model generates paths using myPathGenerator
        // each path is priced using myPathPricer
        // prices will be accumulated into statisticsAccumulator
        MonteCarloModel<SingleVariate,PseudoRandom>
            MCSimulation(myPathGenerator,
                         myPathPricer,
                         statisticsAccumulator,
                         false);

        // the model simulates nSamples paths
        MCSimulation.addSamples(nSamples);

        // the sampleAccumulator method
        // gives access to all the methods of statisticsAccumulator
        Real PLMean  = MCSimulation.sampleAccumulator().mean();
        Real PLStDev = MCSimulation.sampleAccumulator().standardDeviation();
        Real PLSkew  = MCSimulation.sampleAccumulator().skewness();
        Real PLKurt  = MCSimulation.sampleAccumulator().kurtosis();

        // Derman and Kamal's formula
        Real theorStD = std::sqrt(M_PI/4/nTimeSteps)*vega_*sigma_;

        std::cout << std::fixed
                  << std::setw(8) << nSamples << " | "
                  << std::setw(8) << nTimeSteps << " | "
                  << std::setw(8) << std::setprecision(3) << PLMean << " | "
                  << std::setw(8) << std::setprecision(2) << PLStDev << " | "
                  << std::setw(12) << std::setprecision(2) << theorStD << " | "
                  << std::setw(8) << std::setprecision(2) << PLSkew << " | "
                  << std::setw(8) << std::setprecision(2) << PLKurt << std::endl;
    }