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(-r_*maturity_);
84 DiscountFactor qDiscount = 1.0;
85 Real forward = s0_*qDiscount/rDiscount;
86 Real stdDev = std::sqrt(sigma_*sigma_*maturity_);
87 auto payoff = ext::make_shared<PlainVanillaPayoff>(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_);
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(hedgesNum, scenarios);
179
180 hedgesNum = 84;
181 rp.compute(hedgesNum, 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(-r_*maturity_);
225 DiscountFactor qDiscount = std::exp(-stockDividendYield*maturity_);
226 Real forward = stock*qDiscount/rDiscount;
227 Real stdDev = std::sqrt(sigma_*sigma_*maturity_);
228 auto payoff = ext::make_shared<PlainVanillaPayoff>(type_,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(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( 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(-r_*(maturity_-t));
255 qDiscount = std::exp(-stockDividendYield*(maturity_-t));
256 forward = stock*qDiscount/rDiscount;
257 stdDev = std::sqrt(sigma_*sigma_*(maturity_-t));
258 BlackCalculator black(payoff,forward,stdDev,rDiscount);
259
260 // recalculate delta
261 delta = black.delta(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( 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>(s0_));
305 Handle<YieldTermStructure> riskFreeRate(
306 ext::make_shared<FlatForward>(today, r_, dayCount));
307 Handle<YieldTermStructure> dividendYield(
308 ext::make_shared<FlatForward>(today, 0.0, dayCount));
309 Handle<BlackVolTermStructure> volatility(
310 ext::make_shared<BlackConstantVol>(today, calendar, sigma_, dayCount));
311 auto diffusion = ext::make_shared<BlackScholesMertonProcess>(
312 stateVariable, dividendYield, riskFreeRate, 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(nTimeSteps, 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 diffusion, maturity_, nTimeSteps,
325 rsg, 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 payoff_.optionType(), payoff_.strike(),
332 r_, maturity_, 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(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 }