| 1 | /* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ |
| 2 | |
| 3 | /*! |
| 4 | Copyright (C) 2000, 2001, 2002, 2003 RiskMap srl |
| 5 | Copyright (C) 2003, 2004, 2005, 2006, 2007 StatPro Italia srl |
| 6 | |
| 7 | This file is part of QuantLib, a free-software/open-source library |
| 8 | for financial quantitative analysts and developers - http://quantlib.org/ |
| 9 | |
| 10 | QuantLib is free software: you can redistribute it and/or modify it |
| 11 | under the terms of the QuantLib license. You should have received a |
| 12 | copy of the license along with this program; if not, please email |
| 13 | <quantlib-dev@lists.sf.net>. The license is also available online at |
| 14 | <http://quantlib.org/license.shtml>. |
| 15 | |
| 16 | This program is distributed in the hope that it will be useful, but WITHOUT |
| 17 | ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS |
| 18 | FOR A PARTICULAR PURPOSE. See the license for more details. |
| 19 | */ |
| 20 | |
| 21 | /* This example computes profit and loss of a discrete interval hedging |
| 22 | strategy and compares with the results of Derman & Kamal's (Goldman Sachs |
| 23 | Equity Derivatives Research) Research Note: "When You Cannot Hedge |
| 24 | Continuously: The Corrections to Black-Scholes" |
| 25 | http://www.ederman.com/emanuelderman/GSQSpapers/when_you_cannot_hedge.pdf |
| 26 | |
| 27 | Suppose an option hedger sells an European option and receives the |
| 28 | Black-Scholes value as the options premium. |
| 29 | Then he follows a Black-Scholes hedging strategy, rehedging at discrete, |
| 30 | evenly spaced time intervals as the underlying stock changes. At |
| 31 | expiration, the hedger delivers the option payoff to the option holder, |
| 32 | and unwinds the hedge. We are interested in understanding the final |
| 33 | profit or loss of this strategy. |
| 34 | |
| 35 | If the hedger had followed the exact Black-Scholes replication strategy, |
| 36 | re-hedging continuously as the underlying stock evolved towards its final |
| 37 | value at expiration, then, no matter what path the stock took, the final |
| 38 | P&L would be exactly zero. When the replication strategy deviates from |
| 39 | the exact Black-Scholes method, the final P&L may deviate from zero. This |
| 40 | deviation is called the replication error. When the hedger rebalances at |
| 41 | discrete rather than continuous intervals, the hedge is imperfect and the |
| 42 | replication is inexact. The more often hedging occurs, the smaller the |
| 43 | replication error. |
| 44 | |
| 45 | We examine the range of possibilities, computing the replication error. |
| 46 | */ |
| 47 | |
| 48 | #include <ql/qldefines.hpp> |
| 49 | #if !defined(BOOST_ALL_NO_LIB) && defined(BOOST_MSVC) |
| 50 | # include <ql/auto_link.hpp> |
| 51 | #endif |
| 52 | #include <ql/methods/montecarlo/montecarlomodel.hpp> |
| 53 | #include <ql/processes/blackscholesprocess.hpp> |
| 54 | #include <ql/termstructures/yield/flatforward.hpp> |
| 55 | #include <ql/termstructures/volatility/equityfx/blackconstantvol.hpp> |
| 56 | #include <ql/pricingengines/blackcalculator.hpp> |
| 57 | #include <ql/quotes/simplequote.hpp> |
| 58 | #include <ql/time/calendars/target.hpp> |
| 59 | #include <ql/time/daycounters/actual365fixed.hpp> |
| 60 | |
| 61 | #include <iostream> |
| 62 | #include <iomanip> |
| 63 | |
| 64 | using namespace QuantLib; |
| 65 | |
| 66 | /* The ReplicationError class carries out Monte Carlo simulations to evaluate |
| 67 | the outcome (the replication error) of the discrete hedging strategy over |
| 68 | different, randomly generated scenarios of future stock price evolution. |
| 69 | */ |
| 70 | class ReplicationError |
| 71 | { |
| 72 | public: |
| 73 | ReplicationError(Option::Type type, |
| 74 | Time maturity, |
| 75 | Real strike, |
| 76 | Real s0, |
| 77 | Volatility sigma, |
| 78 | Rate r) |
| 79 | : maturity_(maturity), payoff_(type, strike), s0_(s0), |
| 80 | sigma_(sigma), r_(r) { |
| 81 | |
| 82 | // value of the option |
| 83 | DiscountFactor rDiscount = std::exp(x: -r_*maturity_); |
| 84 | DiscountFactor qDiscount = 1.0; |
| 85 | Real forward = s0_*qDiscount/rDiscount; |
| 86 | Real stdDev = std::sqrt(x: sigma_*sigma_*maturity_); |
| 87 | auto payoff = ext::make_shared<PlainVanillaPayoff>(args&: payoff_); |
| 88 | BlackCalculator black(payoff,forward,stdDev,rDiscount); |
| 89 | std::cout << "Option value: " << black.value() << std::endl; |
| 90 | |
| 91 | // store option's vega, since Derman and Kamal's formula needs it |
| 92 | vega_ = black.vega(maturity: maturity_); |
| 93 | |
| 94 | std::cout << std::endl; |
| 95 | |
| 96 | std::cout << std::setw(8) << " " << " | " |
| 97 | << std::setw(8) << " " << " | " |
| 98 | << std::setw(8) << "P&L" << " | " |
| 99 | << std::setw(8) << "P&L" << " | " |
| 100 | << std::setw(12) << "Derman&Kamal" << " | " |
| 101 | << std::setw(8) << "P&L" << " | " |
| 102 | << std::setw(8) << "P&L" << std::endl; |
| 103 | |
| 104 | std::cout << std::setw(8) << "samples" << " | " |
| 105 | << std::setw(8) << "trades" << " | " |
| 106 | << std::setw(8) << "mean" << " | " |
| 107 | << std::setw(8) << "std.dev." << " | " |
| 108 | << std::setw(12) << "formula" << " | " |
| 109 | << std::setw(8) << "skewness" << " | " |
| 110 | << std::setw(8) << "kurtosis" << std::endl; |
| 111 | |
| 112 | std::cout << std::string(78, '-') << std::endl; |
| 113 | } |
| 114 | |
| 115 | // the actual replication error computation |
| 116 | void compute(Size nTimeSteps, Size nSamples); |
| 117 | private: |
| 118 | Time maturity_; |
| 119 | PlainVanillaPayoff payoff_; |
| 120 | Real s0_; |
| 121 | Volatility sigma_; |
| 122 | Rate r_; |
| 123 | Real vega_; |
| 124 | }; |
| 125 | |
| 126 | // The key for the MonteCarlo simulation is to have a PathPricer that |
| 127 | // implements a value(const Path& path) method. |
| 128 | // This method prices the portfolio for each Path of the random variable |
| 129 | class ReplicationPathPricer : public PathPricer<Path> { |
| 130 | public: |
| 131 | // real constructor |
| 132 | ReplicationPathPricer(Option::Type type, |
| 133 | Real strike, |
| 134 | Rate r, |
| 135 | Time maturity, |
| 136 | Volatility sigma) |
| 137 | : type_(type), strike_(strike), |
| 138 | r_(r), maturity_(maturity), sigma_(sigma) { |
| 139 | QL_REQUIRE(strike_ > 0.0, "strike must be positive" ); |
| 140 | QL_REQUIRE(r_ >= 0.0, |
| 141 | "risk free rate (r) must be positive or zero" ); |
| 142 | QL_REQUIRE(maturity_ > 0.0, "maturity must be positive" ); |
| 143 | QL_REQUIRE(sigma_ >= 0.0, |
| 144 | "volatility (sigma) must be positive or zero" ); |
| 145 | |
| 146 | } |
| 147 | // The value() method encapsulates the pricing code |
| 148 | Real operator()(const Path& path) const override; |
| 149 | |
| 150 | private: |
| 151 | Option::Type type_; |
| 152 | Real strike_; |
| 153 | Rate r_; |
| 154 | Time maturity_; |
| 155 | Volatility sigma_; |
| 156 | }; |
| 157 | |
| 158 | |
| 159 | // Compute Replication Error as in the Derman and Kamal's research note |
| 160 | int main(int, char* []) { |
| 161 | |
| 162 | try { |
| 163 | |
| 164 | std::cout << std::endl; |
| 165 | |
| 166 | Time maturity = 1.0/12.0; // 1 month |
| 167 | Real strike = 100; |
| 168 | Real underlying = 100; |
| 169 | Volatility volatility = 0.20; // 20% |
| 170 | Rate riskFreeRate = 0.05; // 5% |
| 171 | ReplicationError rp(Option::Call, maturity, strike, underlying, |
| 172 | volatility, riskFreeRate); |
| 173 | |
| 174 | Size scenarios = 50000; |
| 175 | Size hedgesNum; |
| 176 | |
| 177 | hedgesNum = 21; |
| 178 | rp.compute(nTimeSteps: hedgesNum, nSamples: scenarios); |
| 179 | |
| 180 | hedgesNum = 84; |
| 181 | rp.compute(nTimeSteps: hedgesNum, nSamples: scenarios); |
| 182 | |
| 183 | return 0; |
| 184 | } catch (std::exception& e) { |
| 185 | std::cerr << e.what() << std::endl; |
| 186 | return 1; |
| 187 | } catch (...) { |
| 188 | std::cerr << "unknown error" << std::endl; |
| 189 | return 1; |
| 190 | } |
| 191 | } |
| 192 | |
| 193 | |
| 194 | /* The actual computation of the Profit&Loss for each single path. |
| 195 | |
| 196 | In each scenario N rehedging trades spaced evenly in time over |
| 197 | the life of the option are carried out, using the Black-Scholes |
| 198 | hedge ratio. |
| 199 | */ |
| 200 | Real ReplicationPathPricer::operator()(const Path& path) const { |
| 201 | |
| 202 | Size n = path.length()-1; |
| 203 | QL_REQUIRE(n>0, "the path cannot be empty" ); |
| 204 | |
| 205 | // discrete hedging interval |
| 206 | Time dt = maturity_/n; |
| 207 | |
| 208 | // For simplicity, we assume the stock pays no dividends. |
| 209 | Rate stockDividendYield = 0.0; |
| 210 | |
| 211 | // let's start |
| 212 | Time t = 0; |
| 213 | |
| 214 | // stock value at t=0 |
| 215 | Real stock = path.front(); |
| 216 | |
| 217 | // money account at t=0 |
| 218 | Real money_account = 0.0; |
| 219 | |
| 220 | /************************/ |
| 221 | /*** the initial deal ***/ |
| 222 | /************************/ |
| 223 | // option fair price (Black-Scholes) at t=0 |
| 224 | DiscountFactor rDiscount = std::exp(x: -r_*maturity_); |
| 225 | DiscountFactor qDiscount = std::exp(x: -stockDividendYield*maturity_); |
| 226 | Real forward = stock*qDiscount/rDiscount; |
| 227 | Real stdDev = std::sqrt(x: sigma_*sigma_*maturity_); |
| 228 | auto payoff = ext::make_shared<PlainVanillaPayoff>(args: type_,args: strike_); |
| 229 | BlackCalculator black(payoff,forward,stdDev,rDiscount); |
| 230 | // sell the option, cash in its premium |
| 231 | money_account += black.value(); |
| 232 | // compute delta |
| 233 | Real delta = black.delta(spot: stock); |
| 234 | // delta-hedge the option buying stock |
| 235 | Real stockAmount = delta; |
| 236 | money_account -= stockAmount*stock; |
| 237 | |
| 238 | /**********************************/ |
| 239 | /*** hedging during option life ***/ |
| 240 | /**********************************/ |
| 241 | for (Size step = 0; step < n-1; step++){ |
| 242 | |
| 243 | // time flows |
| 244 | t += dt; |
| 245 | |
| 246 | // accruing on the money account |
| 247 | money_account *= std::exp( x: r_*dt ); |
| 248 | |
| 249 | // stock growth: |
| 250 | stock = path[step+1]; |
| 251 | |
| 252 | // recalculate option value at the current stock value, |
| 253 | // and the current time to maturity |
| 254 | rDiscount = std::exp(x: -r_*(maturity_-t)); |
| 255 | qDiscount = std::exp(x: -stockDividendYield*(maturity_-t)); |
| 256 | forward = stock*qDiscount/rDiscount; |
| 257 | stdDev = std::sqrt(x: sigma_*sigma_*(maturity_-t)); |
| 258 | BlackCalculator black(payoff,forward,stdDev,rDiscount); |
| 259 | |
| 260 | // recalculate delta |
| 261 | delta = black.delta(spot: stock); |
| 262 | |
| 263 | // re-hedging |
| 264 | money_account -= (delta - stockAmount)*stock; |
| 265 | stockAmount = delta; |
| 266 | } |
| 267 | |
| 268 | /*************************/ |
| 269 | /*** option expiration ***/ |
| 270 | /*************************/ |
| 271 | // last accrual on my money account |
| 272 | money_account *= std::exp( x: r_*dt ); |
| 273 | // last stock growth |
| 274 | stock = path[n]; |
| 275 | |
| 276 | // the hedger delivers the option payoff to the option holder |
| 277 | Real optionPayoff = PlainVanillaPayoff(type_, strike_)(stock); |
| 278 | money_account -= optionPayoff; |
| 279 | |
| 280 | // and unwinds the hedge selling his stock position |
| 281 | money_account += stockAmount*stock; |
| 282 | |
| 283 | // final Profit&Loss |
| 284 | return money_account; |
| 285 | } |
| 286 | |
| 287 | |
| 288 | // The computation over nSamples paths of the P&L distribution |
| 289 | void ReplicationError::compute(Size nTimeSteps, Size nSamples) |
| 290 | { |
| 291 | QL_REQUIRE(nTimeSteps>0, "the number of steps must be > 0" ); |
| 292 | |
| 293 | // hedging interval |
| 294 | // Time tau = maturity_ / nTimeSteps; |
| 295 | |
| 296 | /* Black-Scholes framework: the underlying stock price evolves |
| 297 | lognormally with a fixed known volatility that stays constant |
| 298 | throughout time. |
| 299 | */ |
| 300 | Calendar calendar = TARGET(); |
| 301 | Date today = Date::todaysDate(); |
| 302 | DayCounter dayCount = Actual365Fixed(); |
| 303 | Handle<Quote> stateVariable( |
| 304 | ext::make_shared<SimpleQuote>(args&: s0_)); |
| 305 | Handle<YieldTermStructure> riskFreeRate( |
| 306 | ext::make_shared<FlatForward>(args&: today, args&: r_, args&: dayCount)); |
| 307 | Handle<YieldTermStructure> dividendYield( |
| 308 | ext::make_shared<FlatForward>(args&: today, args: 0.0, args&: dayCount)); |
| 309 | Handle<BlackVolTermStructure> volatility( |
| 310 | ext::make_shared<BlackConstantVol>(args&: today, args&: calendar, args&: sigma_, args&: dayCount)); |
| 311 | auto diffusion = ext::make_shared<BlackScholesMertonProcess>( |
| 312 | args&: stateVariable, args&: dividendYield, args&: riskFreeRate, args&: volatility); |
| 313 | |
| 314 | // Black Scholes equation rules the path generator: |
| 315 | // at each step the log of the stock |
| 316 | // will have drift and sigma^2 variance |
| 317 | PseudoRandom::rsg_type rsg = |
| 318 | PseudoRandom::make_sequence_generator(dimension: nTimeSteps, seed: 0); |
| 319 | |
| 320 | bool brownianBridge = false; |
| 321 | |
| 322 | typedef SingleVariate<PseudoRandom>::path_generator_type generator_type; |
| 323 | auto myPathGenerator = ext::make_shared<generator_type>( |
| 324 | args&: diffusion, args&: maturity_, args&: nTimeSteps, |
| 325 | args&: rsg, args&: brownianBridge); |
| 326 | |
| 327 | // The replication strategy's Profit&Loss is computed for each path |
| 328 | // of the stock. The path pricer knows how to price a path using its |
| 329 | // value() method |
| 330 | auto myPathPricer = ext::make_shared<ReplicationPathPricer>( |
| 331 | args: payoff_.optionType(), args: payoff_.strike(), |
| 332 | args&: r_, args&: maturity_, args&: sigma_); |
| 333 | |
| 334 | // a statistics accumulator for the path-dependant Profit&Loss values |
| 335 | Statistics statisticsAccumulator; |
| 336 | |
| 337 | // The Monte Carlo model generates paths using myPathGenerator |
| 338 | // each path is priced using myPathPricer |
| 339 | // prices will be accumulated into statisticsAccumulator |
| 340 | MonteCarloModel<SingleVariate,PseudoRandom> |
| 341 | MCSimulation(myPathGenerator, |
| 342 | myPathPricer, |
| 343 | statisticsAccumulator, |
| 344 | false); |
| 345 | |
| 346 | // the model simulates nSamples paths |
| 347 | MCSimulation.addSamples(samples: nSamples); |
| 348 | |
| 349 | // the sampleAccumulator method |
| 350 | // gives access to all the methods of statisticsAccumulator |
| 351 | Real PLMean = MCSimulation.sampleAccumulator().mean(); |
| 352 | Real PLStDev = MCSimulation.sampleAccumulator().standardDeviation(); |
| 353 | Real PLSkew = MCSimulation.sampleAccumulator().skewness(); |
| 354 | Real PLKurt = MCSimulation.sampleAccumulator().kurtosis(); |
| 355 | |
| 356 | // Derman and Kamal's formula |
| 357 | Real theorStD = std::sqrt(M_PI/4/nTimeSteps)*vega_*sigma_; |
| 358 | |
| 359 | |
| 360 | std::cout << std::fixed |
| 361 | << std::setw(8) << nSamples << " | " |
| 362 | << std::setw(8) << nTimeSteps << " | " |
| 363 | << std::setw(8) << std::setprecision(3) << PLMean << " | " |
| 364 | << std::setw(8) << std::setprecision(2) << PLStDev << " | " |
| 365 | << std::setw(12) << std::setprecision(2) << theorStD << " | " |
| 366 | << std::setw(8) << std::setprecision(2) << PLSkew << " | " |
| 367 | << std::setw(8) << std::setprecision(2) << PLKurt << std::endl; |
| 368 | } |
| 369 | |