DiscreteHedging.cpp

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

00001 /* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
00002 
00021 /*  This example computes profit and loss of a discrete interval hedging
00022     strategy and compares with the results of Derman & Kamal's (Goldman Sachs
00023     Equity Derivatives Research) Research Note: "When You Cannot Hedge
00024     Continuously: The Corrections to Black-Scholes"
00025     http://www.ederman.com/emanuelderman/GSQSpapers/when_you_cannot_hedge.pdf
00026 
00027     Suppose an option hedger sells an European option and receives the
00028     Black-Scholes value as the options premium.
00029     Then he follows a Black-Scholes hedging strategy, rehedging at discrete,
00030     evenly spaced time intervals as the underlying stock changes. At
00031     expiration, the hedger delivers the option payoff to the option holder,
00032     and unwinds the hedge. We are interested in understanding the final
00033     profit or loss of this strategy.
00034 
00035     If the hedger had followed the exact Black-Scholes replication strategy,
00036     re-hedging continuously as the underlying stock evolved towards its final
00037     value at expiration, then, no matter what path the stock took, the final
00038     P&L would be exactly zero. When the replication strategy deviates from
00039     the exact Black-Scholes method, the final P&L may deviate from zero. This
00040     deviation is called the replication error. When the hedger rebalances at
00041     discrete rather than continuous intervals, the hedge is imperfect and the
00042     replication is inexact. The more often hedging occurs, the smaller the
00043     replication error.
00044 
00045     We examine the range of possibilities, computing the replication error.
00046 */
00047 
00048 // the only header you need to use QuantLib
00049 #define BOOST_LIB_DIAGNOSTIC
00050 #  include <ql/quantlib.hpp>
00051 #undef BOOST_LIB_DIAGNOSTIC
00052 
00053 #ifdef BOOST_MSVC
00054 /* Uncomment the following lines to unmask floating-point
00055    exceptions. Warning: unpredictable results can arise...
00056 
00057    See http://www.wilmott.com/messageview.cfm?catid=10&threadid=9481
00058    Is there anyone with a definitive word about this?
00059 */
00060 // #include <float.h>
00061 // namespace { unsigned int u = _controlfp(_EM_INEXACT, _MCW_EM); }
00062 #endif
00063 
00064 #include <boost/timer.hpp>
00065 #include <iostream>
00066 #include <iomanip>
00067 
00068 using namespace QuantLib;
00069 
00070 #if defined(QL_ENABLE_SESSIONS)
00071 namespace QuantLib {
00072 
00073     Integer sessionId() { return 0; }
00074 
00075 }
00076 #endif
00077 
00078 
00079 /* The ReplicationError class carries out Monte Carlo simulations to evaluate
00080    the outcome (the replication error) of the discrete hedging strategy over
00081    different, randomly generated scenarios of future stock price evolution.
00082 */
00083 class ReplicationError
00084 {
00085 public:
00086     ReplicationError(Option::Type type,
00087                      Time maturity,
00088                      Real strike,
00089                      Real s0,
00090                      Volatility sigma,
00091                      Rate r)
00092     : maturity_(maturity), payoff_(type, strike), s0_(s0),
00093       sigma_(sigma), r_(r) {
00094 
00095         // value of the option
00096         DiscountFactor rDiscount = std::exp(-r_*maturity_);
00097         DiscountFactor qDiscount = 1.0;
00098         Real forward = s0_*qDiscount/rDiscount;
00099         Real stdDev = std::sqrt(sigma_*sigma_*maturity_);
00100         boost::shared_ptr<StrikedTypePayoff> payoff(
00101                                              new PlainVanillaPayoff(payoff_));
00102         BlackCalculator black(payoff,forward,stdDev,rDiscount);
00103         std::cout << "Option value: " << black.value() << std::endl;
00104 
00105         // store option's vega, since Derman and Kamal's formula needs it
00106         vega_ = black.vega(maturity_);
00107 
00108         std::cout << std::endl;
00109         std::cout <<
00110             "        |        | P&L  \t|  P&L    | Derman&Kamal | P&L"
00111             "      \t| P&L" << std::endl;
00112 
00113         std::cout <<
00114             "samples | trades | Mean \t| Std Dev | Formula      |"
00115             " skewness \t| kurt." << std::endl;
00116 
00117         std::cout << "---------------------------------"
00118             "----------------------------------------------" << std::endl;
00119     }
00120 
00121     // the actual replication error computation
00122     void compute(Size nTimeSteps, Size nSamples);
00123 private:
00124     Time maturity_;
00125     PlainVanillaPayoff payoff_;
00126     Real s0_;
00127     Volatility sigma_;
00128     Rate r_;
00129     Real vega_;
00130 };
00131 
00132 // The key for the MonteCarlo simulation is to have a PathPricer that
00133 // implements a value(const Path& path) method.
00134 // This method prices the portfolio for each Path of the random variable
00135 class ReplicationPathPricer : public PathPricer<Path> {
00136   public:
00137     // real constructor
00138     ReplicationPathPricer(Option::Type type,
00139                           Real strike,
00140                           Rate r,
00141                           Time maturity,
00142                           Volatility sigma)
00143     : type_(type), strike_(strike),
00144       r_(r), maturity_(maturity), sigma_(sigma) {
00145         QL_REQUIRE(strike_ > 0.0, "strike must be positive");
00146         QL_REQUIRE(r_ >= 0.0,
00147                    "risk free rate (r) must be positive or zero");
00148         QL_REQUIRE(maturity_ > 0.0, "maturity must be positive");
00149         QL_REQUIRE(sigma_ >= 0.0,
00150                    "volatility (sigma) must be positive or zero");
00151 
00152     }
00153     // The value() method encapsulates the pricing code
00154     Real operator()(const Path& path) const;
00155 
00156   private:
00157     Option::Type type_;
00158     Real strike_;
00159     Rate r_;
00160     Time maturity_;
00161     Volatility sigma_;
00162 };
00163 
00164 
00165 // Compute Replication Error as in the Derman and Kamal's research note
00166 int main(int, char* []) {
00167 
00168     try {
00169 
00170         boost::timer timer;
00171         std::cout << std::endl;
00172 
00173         Time maturity = 1.0/12.0;   // 1 month
00174         Real strike = 100;
00175         Real underlying = 100;
00176         Volatility volatility = 0.20; // 20%
00177         Rate riskFreeRate = 0.05; // 5%
00178         ReplicationError rp(Option::Call, maturity, strike, underlying,
00179                 volatility, riskFreeRate);
00180 
00181         Size scenarios = 50000;
00182         Size hedgesNum;
00183 
00184         hedgesNum = 21;
00185         rp.compute(hedgesNum, scenarios);
00186 
00187         hedgesNum = 84;
00188         rp.compute(hedgesNum, scenarios);
00189 
00190         Real seconds = timer.elapsed();
00191         Integer hours = int(seconds/3600);
00192         seconds -= hours * 3600;
00193         Integer minutes = int(seconds/60);
00194         seconds -= minutes * 60;
00195         std::cout << " \nRun completed in ";
00196         if (hours > 0)
00197             std::cout << hours << " h ";
00198         if (hours > 0 || minutes > 0)
00199             std::cout << minutes << " m ";
00200         std::cout << std::fixed << std::setprecision(0)
00201                   << seconds << " s\n" << std::endl;
00202 
00203         return 0;
00204     } catch (std::exception& e) {
00205         std::cout << e.what() << std::endl;
00206         return 1;
00207     } catch (...) {
00208         std::cout << "unknown error" << std::endl;
00209         return 1;
00210     }
00211 }
00212 
00213 
00214 /* The actual computation of the Profit&Loss for each single path.
00215 
00216    In each scenario N rehedging trades spaced evenly in time over
00217    the life of the option are carried out, using the Black-Scholes
00218    hedge ratio.
00219 */
00220 Real ReplicationPathPricer::operator()(const Path& path) const {
00221 
00222     Size n = path.length()-1;
00223     QL_REQUIRE(n>0, "the path cannot be empty");
00224 
00225     // discrete hedging interval
00226     Time dt = maturity_/n;
00227 
00228     // For simplicity, we assume the stock pays no dividends.
00229     Rate stockDividendYield = 0.0;
00230 
00231     // let's start
00232     Time t = 0;
00233 
00234     // stock value at t=0
00235     Real stock = path.front();
00236 
00237     // money account at t=0
00238     Real money_account = 0.0;
00239 
00240     /************************/
00241     /*** the initial deal ***/
00242     /************************/
00243     // option fair price (Black-Scholes) at t=0
00244     DiscountFactor rDiscount = std::exp(-r_*maturity_);
00245     DiscountFactor qDiscount = std::exp(-stockDividendYield*maturity_);
00246     Real forward = stock*qDiscount/rDiscount;
00247     Real stdDev = std::sqrt(sigma_*sigma_*maturity_);
00248     boost::shared_ptr<StrikedTypePayoff> payoff(
00249                                        new PlainVanillaPayoff(type_,strike_));
00250     BlackCalculator black(payoff,forward,stdDev,rDiscount);
00251     // sell the option, cash in its premium
00252     money_account += black.value();
00253     // compute delta
00254     Real delta = black.delta(stock);
00255     // delta-hedge the option buying stock
00256     Real stockAmount = delta;
00257     money_account -= stockAmount*stock;
00258 
00259     /**********************************/
00260     /*** hedging during option life ***/
00261     /**********************************/
00262     for (Size step = 0; step < n-1; step++){
00263 
00264         // time flows
00265         t += dt;
00266 
00267         // accruing on the money account
00268         money_account *= std::exp( r_*dt );
00269 
00270         // stock growth:
00271         stock = path[step+1];
00272 
00273         // recalculate option value at the current stock value,
00274         // and the current time to maturity
00275         rDiscount = std::exp(-r_*(maturity_-t));
00276         qDiscount = std::exp(-stockDividendYield*(maturity_-t));
00277         forward = stock*qDiscount/rDiscount;
00278         stdDev = std::sqrt(sigma_*sigma_*(maturity_-t));
00279         BlackCalculator black(payoff,forward,stdDev,rDiscount);
00280 
00281         // recalculate delta
00282         delta = black.delta(stock);
00283 
00284         // re-hedging
00285         money_account -= (delta - stockAmount)*stock;
00286         stockAmount = delta;
00287     }
00288 
00289     /*************************/
00290     /*** option expiration ***/
00291     /*************************/
00292     // last accrual on my money account
00293     money_account *= std::exp( r_*dt );
00294     // last stock growth
00295     stock = path[n];
00296 
00297     // the hedger delivers the option payoff to the option holder
00298     Real optionPayoff = PlainVanillaPayoff(type_, strike_)(stock);
00299     money_account -= optionPayoff;
00300 
00301     // and unwinds the hedge selling his stock position
00302     money_account += stockAmount*stock;
00303 
00304     // final Profit&Loss
00305     return money_account;
00306 }
00307 
00308 
00309 // The computation over nSamples paths of the P&L distribution
00310 void ReplicationError::compute(Size nTimeSteps, Size nSamples)
00311 {
00312     QL_REQUIRE(nTimeSteps>0, "the number of steps must be > 0");
00313 
00314     // hedging interval
00315     // Time tau = maturity_ / nTimeSteps;
00316 
00317     /* Black-Scholes framework: the underlying stock price evolves
00318        lognormally with a fixed known volatility that stays constant
00319        throughout time.
00320     */
00321     Date today = Date::todaysDate();
00322     DayCounter dayCount = Actual365Fixed();
00323     Handle<Quote> stateVariable(
00324                           boost::shared_ptr<Quote>(new SimpleQuote(s0_)));
00325     Handle<YieldTermStructure> riskFreeRate(
00326                           boost::shared_ptr<YieldTermStructure>(
00327                                       new FlatForward(today, r_, dayCount)));
00328     Handle<YieldTermStructure> dividendYield(
00329                           boost::shared_ptr<YieldTermStructure>(
00330                                       new FlatForward(today, 0.0, dayCount)));
00331     Handle<BlackVolTermStructure> volatility(
00332                           boost::shared_ptr<BlackVolTermStructure>(
00333                                new BlackConstantVol(today, sigma_, dayCount)));
00334     boost::shared_ptr<StochasticProcess1D> diffusion(
00335                    new BlackScholesMertonProcess(stateVariable, dividendYield,
00336                                                  riskFreeRate, volatility));
00337 
00338     // Black Scholes equation rules the path generator:
00339     // at each step the log of the stock
00340     // will have drift and sigma^2 variance
00341     PseudoRandom::rsg_type rsg =
00342         PseudoRandom::make_sequence_generator(nTimeSteps, 0);
00343 
00344     bool brownianBridge = false;
00345 
00346     typedef SingleVariate<PseudoRandom>::path_generator_type generator_type;
00347     boost::shared_ptr<generator_type> myPathGenerator(new
00348         generator_type(diffusion, maturity_, nTimeSteps,
00349                        rsg, brownianBridge));
00350 
00351     // The replication strategy's Profit&Loss is computed for each path
00352     // of the stock. The path pricer knows how to price a path using its
00353     // value() method
00354     boost::shared_ptr<PathPricer<Path> > myPathPricer(new
00355         ReplicationPathPricer(payoff_.optionType(), payoff_.strike(),
00356                               r_, maturity_, sigma_));
00357 
00358     // a statistics accumulator for the path-dependant Profit&Loss values
00359     Statistics statisticsAccumulator;
00360 
00361     // The Monte Carlo model generates paths using myPathGenerator
00362     // each path is priced using myPathPricer
00363     // prices will be accumulated into statisticsAccumulator
00364     MonteCarloModel<SingleVariate,PseudoRandom>
00365         MCSimulation(myPathGenerator,
00366                      myPathPricer,
00367                      statisticsAccumulator,
00368                      false);
00369 
00370     // the model simulates nSamples paths
00371     MCSimulation.addSamples(nSamples);
00372 
00373     // the sampleAccumulator method
00374     // gives access to all the methods of statisticsAccumulator
00375     Real PLMean  = MCSimulation.sampleAccumulator().mean();
00376     Real PLStDev = MCSimulation.sampleAccumulator().standardDeviation();
00377     Real PLSkew  = MCSimulation.sampleAccumulator().skewness();
00378     Real PLKurt  = MCSimulation.sampleAccumulator().kurtosis();
00379 
00380     // Derman and Kamal's formula
00381     Real theorStD = std::sqrt(M_PI/4/nTimeSteps)*vega_*sigma_;
00382 
00383 
00384     std::cout << std::fixed
00385               << nSamples << "\t| "
00386               << nTimeSteps << "\t | "
00387               << std::setprecision(3) << PLMean << " \t| "
00388               << std::setprecision(2) << PLStDev << " \t  | "
00389               << std::setprecision(2) << theorStD << " \t | "
00390               << std::setprecision(2) << PLSkew << " \t| "
00391               << std::setprecision(2) << PLKurt << std::endl;
00392 }