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 }