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 }