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