1 /* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2 
3 /*
4  Copyright (C) 2003, 2007 Ferdinando Ametrano
5  Copyright (C) 2003, 2007 StatPro Italia srl
6  Copyright (C) 2009 Klaus Spanderen
7 
8  This file is part of QuantLib, a free-software/open-source library
9  for financial quantitative analysts and developers - http://quantlib.org/
10 
11  QuantLib is free software: you can redistribute it and/or modify it
12  under the terms of the QuantLib license.  You should have received a
13  copy of the license along with this program; if not, please email
14  <quantlib-dev@lists.sf.net>. The license is also available online at
15  <http://quantlib.org/license.shtml>.
16 
17  This program is distributed in the hope that it will be useful, but WITHOUT
18  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
19  FOR A PARTICULAR PURPOSE.  See the license for more details.
20 */
21 
22 #include "europeanoption.hpp"
23 #include "utilities.hpp"
24 #include <ql/time/calendars/target.hpp>
25 #include <ql/time/daycounters/actual360.hpp>
26 #include <ql/instruments/europeanoption.hpp>
27 #include <ql/math/randomnumbers/rngtraits.hpp>
28 #include <ql/math/interpolations/bicubicsplineinterpolation.hpp>
29 #include <ql/math/interpolations/bilinearinterpolation.hpp>
30 #include <ql/pricingengines/vanilla/analyticeuropeanengine.hpp>
31 #include <ql/pricingengines/vanilla/binomialengine.hpp>
32 #include <ql/pricingengines/vanilla/fdblackscholesvanillaengine.hpp>
33 #include <ql/experimental/variancegamma/fftvanillaengine.hpp>
34 #include <ql/pricingengines/vanilla/mceuropeanengine.hpp>
35 #include <ql/pricingengines/vanilla/integralengine.hpp>
36 #include <ql/termstructures/yield/flatforward.hpp>
37 #include <ql/termstructures/yield/zerocurve.hpp>
38 #include <ql/termstructures/yield/forwardcurve.hpp>
39 #include <ql/termstructures/volatility/equityfx/blackconstantvol.hpp>
40 #include <ql/termstructures/volatility/equityfx/blackvariancesurface.hpp>
41 #include <ql/utilities/dataformatters.hpp>
42 #include <map>
43 
44 using namespace QuantLib;
45 using namespace boost::unit_test_framework;
46 
47 #undef REPORT_FAILURE
48 #define REPORT_FAILURE(greekName, payoff, exercise, s, q, r, today, \
49                        v, expected, calculated, error, tolerance) \
50     BOOST_ERROR(exerciseTypeToString(exercise) << " " \
51                << payoff->optionType() << " option with " \
52                << payoffTypeToString(payoff) << " payoff:\n" \
53                << "    spot value:       " << s << "\n" \
54                << "    strike:           " << payoff->strike() << "\n" \
55                << "    dividend yield:   " << io::rate(q) << "\n" \
56                << "    risk-free rate:   " << io::rate(r) << "\n" \
57                << "    reference date:   " << today << "\n" \
58                << "    maturity:         " << exercise->lastDate() << "\n" \
59                << "    volatility:       " << io::volatility(v) << "\n\n" \
60                << "    expected " << greekName << ":   " << expected << "\n" \
61                << "    calculated " << greekName << ": " << calculated << "\n"\
62                << "    error:            " << error << "\n" \
63                << "    tolerance:        " << tolerance);
64 
65 namespace european_option_test {
66 
67     // utilities
68 
69     struct EuropeanOptionData {
70         Option::Type type;
71         Real strike;
72         Real s;        // spot
73         Rate q;        // dividend
74         Rate r;        // risk-free rate
75         Time t;        // time to maturity
76         Volatility v;  // volatility
77         Real result;   // expected result
78         Real tol;      // tolerance
79     };
80 
81     enum EngineType { Analytic,
82                       JR, CRR, EQP, TGEO, TIAN, LR, JOSHI,
83                       FiniteDifferences,
84                       Integral,
85                       PseudoMonteCarlo, QuasiMonteCarlo,
86                       FFT };
87 
88     ext::shared_ptr<GeneralizedBlackScholesProcess>
makeProcess(const ext::shared_ptr<Quote> & u,const ext::shared_ptr<YieldTermStructure> & q,const ext::shared_ptr<YieldTermStructure> & r,const ext::shared_ptr<BlackVolTermStructure> & vol)89     makeProcess(const ext::shared_ptr<Quote>& u,
90                 const ext::shared_ptr<YieldTermStructure>& q,
91                 const ext::shared_ptr<YieldTermStructure>& r,
92                 const ext::shared_ptr<BlackVolTermStructure>& vol) {
93         return ext::make_shared<BlackScholesMertonProcess>(
94            Handle<Quote>(u),
95                                          Handle<YieldTermStructure>(q),
96                                          Handle<YieldTermStructure>(r),
97                                          Handle<BlackVolTermStructure>(vol));
98     }
99 
100     ext::shared_ptr<VanillaOption>
makeOption(const ext::shared_ptr<StrikedTypePayoff> & payoff,const ext::shared_ptr<Exercise> & exercise,const ext::shared_ptr<Quote> & u,const ext::shared_ptr<YieldTermStructure> & q,const ext::shared_ptr<YieldTermStructure> & r,const ext::shared_ptr<BlackVolTermStructure> & vol,EngineType engineType,Size binomialSteps,Size samples)101     makeOption(const ext::shared_ptr<StrikedTypePayoff>& payoff,
102                const ext::shared_ptr<Exercise>& exercise,
103                const ext::shared_ptr<Quote>& u,
104                const ext::shared_ptr<YieldTermStructure>& q,
105                const ext::shared_ptr<YieldTermStructure>& r,
106                const ext::shared_ptr<BlackVolTermStructure>& vol,
107                EngineType engineType,
108                Size binomialSteps,
109                Size samples) {
110 
111         ext::shared_ptr<GeneralizedBlackScholesProcess> stochProcess =
112             makeProcess(u,q,r,vol);
113 
114         ext::shared_ptr<PricingEngine> engine;
115         switch (engineType) {
116           case Analytic:
117             engine = ext::shared_ptr<PricingEngine>(
118                                     new AnalyticEuropeanEngine(stochProcess));
119             break;
120           case JR:
121             engine = ext::shared_ptr<PricingEngine>(
122                         new BinomialVanillaEngine<JarrowRudd>(stochProcess,
123                                                               binomialSteps));
124             break;
125           case CRR:
126             engine = ext::shared_ptr<PricingEngine>(
127                  new BinomialVanillaEngine<CoxRossRubinstein>(stochProcess,
128                                                               binomialSteps));
129             break;
130           case EQP:
131             engine = ext::shared_ptr<PricingEngine>(
132                  new BinomialVanillaEngine<AdditiveEQPBinomialTree>(
133                                                               stochProcess,
134                                                               binomialSteps));
135             break;
136           case TGEO:
137             engine = ext::shared_ptr<PricingEngine>(
138                         new BinomialVanillaEngine<Trigeorgis>(stochProcess,
139                                                               binomialSteps));
140             break;
141           case TIAN:
142             engine = ext::shared_ptr<PricingEngine>(
143                 new BinomialVanillaEngine<Tian>(stochProcess, binomialSteps));
144             break;
145           case LR:
146             engine = ext::shared_ptr<PricingEngine>(
147                       new BinomialVanillaEngine<LeisenReimer>(stochProcess,
148                                                               binomialSteps));
149             break;
150           case JOSHI:
151             engine = ext::shared_ptr<PricingEngine>(
152               new BinomialVanillaEngine<Joshi4>(stochProcess, binomialSteps));
153             break;
154           case FiniteDifferences:
155             engine = ext::shared_ptr<PricingEngine>(
156                             new FdBlackScholesVanillaEngine(stochProcess,
157                                                             binomialSteps,
158                                                             samples));
159             break;
160           case Integral:
161             engine = ext::shared_ptr<PricingEngine>(
162                                             new IntegralEngine(stochProcess));
163             break;
164           case PseudoMonteCarlo:
165             engine = MakeMCEuropeanEngine<PseudoRandom>(stochProcess)
166                 .withSteps(1)
167                 .withSamples(samples)
168                 .withSeed(42);
169             break;
170           case QuasiMonteCarlo:
171             engine = MakeMCEuropeanEngine<LowDiscrepancy>(stochProcess)
172                 .withSteps(1)
173                 .withSamples(samples);
174             break;
175           case FFT:
176               engine = ext::shared_ptr<PricingEngine>(
177                                           new FFTVanillaEngine(stochProcess));
178             break;
179           default:
180             QL_FAIL("unknown engine type");
181         }
182 
183         ext::shared_ptr<VanillaOption> option(
184             new EuropeanOption(payoff, exercise));
185         option->setPricingEngine(engine);
186         return option;
187     }
188 
timeToDays(Time t)189     Integer timeToDays(Time t) {
190         // FLOATING_POINT_EXCEPTION
191         return Integer(t*360+0.5);
192     }
193 
194 }
195 
196 
testValues()197 void EuropeanOptionTest::testValues() {
198 
199     BOOST_TEST_MESSAGE("Testing European option values...");
200 
201     using namespace european_option_test;
202 
203     SavedSettings backup;
204 
205     /* The data below are from
206        "Option pricing formulas", E.G. Haug, McGraw-Hill 1998
207     */
208     EuropeanOptionData values[] = {
209       // pag 2-8
210       //        type, strike,   spot,    q,    r,    t,  vol,   value,    tol
211       { Option::Call,  65.00,  60.00, 0.00, 0.08, 0.25, 0.30,  2.1334, 1.0e-4},
212       { Option::Put,   95.00, 100.00, 0.05, 0.10, 0.50, 0.20,  2.4648, 1.0e-4},
213       { Option::Put,   19.00,  19.00, 0.10, 0.10, 0.75, 0.28,  1.7011, 1.0e-4},
214       { Option::Call,  19.00,  19.00, 0.10, 0.10, 0.75, 0.28,  1.7011, 1.0e-4},
215       { Option::Call,   1.60,   1.56, 0.08, 0.06, 0.50, 0.12,  0.0291, 1.0e-4},
216       { Option::Put,   70.00,  75.00, 0.05, 0.10, 0.50, 0.35,  4.0870, 1.0e-4},
217       // pag 24
218       { Option::Call, 100.00,  90.00, 0.10, 0.10, 0.10, 0.15,  0.0205, 1.0e-4},
219       { Option::Call, 100.00, 100.00, 0.10, 0.10, 0.10, 0.15,  1.8734, 1.0e-4},
220       { Option::Call, 100.00, 110.00, 0.10, 0.10, 0.10, 0.15,  9.9413, 1.0e-4},
221       { Option::Call, 100.00,  90.00, 0.10, 0.10, 0.10, 0.25,  0.3150, 1.0e-4},
222       { Option::Call, 100.00, 100.00, 0.10, 0.10, 0.10, 0.25,  3.1217, 1.0e-4},
223       { Option::Call, 100.00, 110.00, 0.10, 0.10, 0.10, 0.25, 10.3556, 1.0e-4},
224       { Option::Call, 100.00,  90.00, 0.10, 0.10, 0.10, 0.35,  0.9474, 1.0e-4},
225       { Option::Call, 100.00, 100.00, 0.10, 0.10, 0.10, 0.35,  4.3693, 1.0e-4},
226       { Option::Call, 100.00, 110.00, 0.10, 0.10, 0.10, 0.35, 11.1381, 1.0e-4},
227       { Option::Call, 100.00,  90.00, 0.10, 0.10, 0.50, 0.15,  0.8069, 1.0e-4},
228       { Option::Call, 100.00, 100.00, 0.10, 0.10, 0.50, 0.15,  4.0232, 1.0e-4},
229       { Option::Call, 100.00, 110.00, 0.10, 0.10, 0.50, 0.15, 10.5769, 1.0e-4},
230       { Option::Call, 100.00,  90.00, 0.10, 0.10, 0.50, 0.25,  2.7026, 1.0e-4},
231       { Option::Call, 100.00, 100.00, 0.10, 0.10, 0.50, 0.25,  6.6997, 1.0e-4},
232       { Option::Call, 100.00, 110.00, 0.10, 0.10, 0.50, 0.25, 12.7857, 1.0e-4},
233       { Option::Call, 100.00,  90.00, 0.10, 0.10, 0.50, 0.35,  4.9329, 1.0e-4},
234       { Option::Call, 100.00, 100.00, 0.10, 0.10, 0.50, 0.35,  9.3679, 1.0e-4},
235       { Option::Call, 100.00, 110.00, 0.10, 0.10, 0.50, 0.35, 15.3086, 1.0e-4},
236       { Option::Put,  100.00,  90.00, 0.10, 0.10, 0.10, 0.15,  9.9210, 1.0e-4},
237       { Option::Put,  100.00, 100.00, 0.10, 0.10, 0.10, 0.15,  1.8734, 1.0e-4},
238       { Option::Put,  100.00, 110.00, 0.10, 0.10, 0.10, 0.15,  0.0408, 1.0e-4},
239       { Option::Put,  100.00,  90.00, 0.10, 0.10, 0.10, 0.25, 10.2155, 1.0e-4},
240       { Option::Put,  100.00, 100.00, 0.10, 0.10, 0.10, 0.25,  3.1217, 1.0e-4},
241       { Option::Put,  100.00, 110.00, 0.10, 0.10, 0.10, 0.25,  0.4551, 1.0e-4},
242       { Option::Put,  100.00,  90.00, 0.10, 0.10, 0.10, 0.35, 10.8479, 1.0e-4},
243       { Option::Put,  100.00, 100.00, 0.10, 0.10, 0.10, 0.35,  4.3693, 1.0e-4},
244       { Option::Put,  100.00, 110.00, 0.10, 0.10, 0.10, 0.35,  1.2376, 1.0e-4},
245       { Option::Put,  100.00,  90.00, 0.10, 0.10, 0.50, 0.15, 10.3192, 1.0e-4},
246       { Option::Put,  100.00, 100.00, 0.10, 0.10, 0.50, 0.15,  4.0232, 1.0e-4},
247       { Option::Put,  100.00, 110.00, 0.10, 0.10, 0.50, 0.15,  1.0646, 1.0e-4},
248       { Option::Put,  100.00,  90.00, 0.10, 0.10, 0.50, 0.25, 12.2149, 1.0e-4},
249       { Option::Put,  100.00, 100.00, 0.10, 0.10, 0.50, 0.25,  6.6997, 1.0e-4},
250       { Option::Put,  100.00, 110.00, 0.10, 0.10, 0.50, 0.25,  3.2734, 1.0e-4},
251       { Option::Put,  100.00,  90.00, 0.10, 0.10, 0.50, 0.35, 14.4452, 1.0e-4},
252       { Option::Put,  100.00, 100.00, 0.10, 0.10, 0.50, 0.35,  9.3679, 1.0e-4},
253       { Option::Put,  100.00, 110.00, 0.10, 0.10, 0.50, 0.35,  5.7963, 1.0e-4},
254       // pag 27
255       { Option::Call,  40.00,  42.00, 0.08, 0.04, 0.75, 0.35,  5.0975, 1.0e-4}
256     };
257 
258     DayCounter dc = Actual360();
259     Date today = Date::todaysDate();
260 
261     ext::shared_ptr<SimpleQuote> spot(new SimpleQuote(0.0));
262     ext::shared_ptr<SimpleQuote> qRate(new SimpleQuote(0.0));
263     ext::shared_ptr<YieldTermStructure> qTS = flatRate(today, qRate, dc);
264     ext::shared_ptr<SimpleQuote> rRate(new SimpleQuote(0.0));
265     ext::shared_ptr<YieldTermStructure> rTS = flatRate(today, rRate, dc);
266     ext::shared_ptr<SimpleQuote> vol(new SimpleQuote(0.0));
267     ext::shared_ptr<BlackVolTermStructure> volTS = flatVol(today, vol, dc);
268 
269     for (Size i=0; i<LENGTH(values); i++) {
270 
271         ext::shared_ptr<StrikedTypePayoff> payoff(new
272             PlainVanillaPayoff(values[i].type, values[i].strike));
273         Date exDate = today + timeToDays(values[i].t);
274         ext::shared_ptr<Exercise> exercise(new EuropeanExercise(exDate));
275 
276         spot ->setValue(values[i].s);
277         qRate->setValue(values[i].q);
278         rRate->setValue(values[i].r);
279         vol  ->setValue(values[i].v);
280 
281         ext::shared_ptr<BlackScholesMertonProcess> stochProcess(new
282             BlackScholesMertonProcess(Handle<Quote>(spot),
283                                       Handle<YieldTermStructure>(qTS),
284                                       Handle<YieldTermStructure>(rTS),
285                                       Handle<BlackVolTermStructure>(volTS)));
286         ext::shared_ptr<PricingEngine> engine(
287                                     new AnalyticEuropeanEngine(stochProcess));
288 
289         EuropeanOption option(payoff, exercise);
290         option.setPricingEngine(engine);
291 
292         Real calculated = option.NPV();
293         Real error = std::fabs(calculated-values[i].result);
294         Real tolerance = values[i].tol;
295         if (error>tolerance) {
296             REPORT_FAILURE("value", payoff, exercise, values[i].s,
297                            values[i].q, values[i].r, today,
298                            values[i].v, values[i].result, calculated,
299                            error, tolerance);
300         }
301 
302         engine = ext::shared_ptr<PricingEngine>(
303                     new FdBlackScholesVanillaEngine(stochProcess,200,400));
304         option.setPricingEngine(engine);
305 
306         calculated = option.NPV();
307         error = std::fabs(calculated-values[i].result);
308         tolerance = 1.0e-3;
309         if (error>tolerance) {
310             REPORT_FAILURE("value", payoff, exercise, values[i].s,
311                            values[i].q, values[i].r, today,
312                            values[i].v, values[i].result, calculated,
313                            error, tolerance);
314         }
315     }
316 
317 }
318 
319 
320 
testGreekValues()321 void EuropeanOptionTest::testGreekValues() {
322 
323     BOOST_TEST_MESSAGE("Testing European option greek values...");
324 
325     using namespace european_option_test;
326 
327     SavedSettings backup;
328 
329     /* The data below are from
330        "Option pricing formulas", E.G. Haug, McGraw-Hill 1998
331        pag 11-16
332     */
333     EuropeanOptionData values[] = {
334       //        type, strike,   spot,    q,    r,        t,  vol,  value
335       // delta
336       { Option::Call, 100.00, 105.00, 0.10, 0.10, 0.500000, 0.36,  0.5946, 0 },
337       { Option::Put,  100.00, 105.00, 0.10, 0.10, 0.500000, 0.36, -0.3566, 0 },
338       // elasticity
339       { Option::Put,  100.00, 105.00, 0.10, 0.10, 0.500000, 0.36, -4.8775, 0 },
340       // gamma
341       { Option::Call,  60.00,  55.00, 0.00, 0.10, 0.750000, 0.30,  0.0278, 0 },
342       { Option::Put,   60.00,  55.00, 0.00, 0.10, 0.750000, 0.30,  0.0278, 0 },
343       // vega
344       { Option::Call,  60.00,  55.00, 0.00, 0.10, 0.750000, 0.30, 18.9358, 0 },
345       { Option::Put,   60.00,  55.00, 0.00, 0.10, 0.750000, 0.30, 18.9358, 0 },
346       // theta
347       { Option::Put,  405.00, 430.00, 0.05, 0.07, 1.0/12.0, 0.20,-31.1924, 0 },
348       // theta per day
349       { Option::Put,  405.00, 430.00, 0.05, 0.07, 1.0/12.0, 0.20, -0.0855, 0 },
350       // rho
351       { Option::Call,  75.00,  72.00, 0.00, 0.09, 1.000000, 0.19, 38.7325, 0 },
352       // dividendRho
353       { Option::Put,  490.00, 500.00, 0.05, 0.08, 0.250000, 0.15, 42.2254, 0 }
354     };
355 
356     DayCounter dc = Actual360();
357     Date today = Date::todaysDate();
358 
359     ext::shared_ptr<SimpleQuote> spot(new SimpleQuote(0.0));
360     ext::shared_ptr<SimpleQuote> qRate(new SimpleQuote(0.0));
361     ext::shared_ptr<YieldTermStructure> qTS = flatRate(today, qRate, dc);
362     ext::shared_ptr<SimpleQuote> rRate(new SimpleQuote(0.0));
363     ext::shared_ptr<YieldTermStructure> rTS = flatRate(today, rRate, dc);
364     ext::shared_ptr<SimpleQuote> vol(new SimpleQuote(0.0));
365     ext::shared_ptr<BlackVolTermStructure> volTS = flatVol(today, vol, dc);
366     ext::shared_ptr<BlackScholesMertonProcess> stochProcess(new
367         BlackScholesMertonProcess(Handle<Quote>(spot),
368                                   Handle<YieldTermStructure>(qTS),
369                                   Handle<YieldTermStructure>(rTS),
370                                   Handle<BlackVolTermStructure>(volTS)));
371     ext::shared_ptr<PricingEngine> engine(
372                                     new AnalyticEuropeanEngine(stochProcess));
373 
374     ext::shared_ptr<StrikedTypePayoff> payoff;
375     Date exDate;
376     ext::shared_ptr<Exercise> exercise;
377     ext::shared_ptr<VanillaOption> option;
378     Real calculated;
379 
380     Integer i = -1;
381 
382     i++;
383     payoff = ext::shared_ptr<StrikedTypePayoff>(new
384         PlainVanillaPayoff(values[i].type, values[i].strike));
385     exDate = today + timeToDays(values[i].t);
386     exercise = ext::shared_ptr<Exercise>(new EuropeanExercise(exDate));
387     spot ->setValue(values[i].s);
388     qRate->setValue(values[i].q);
389     rRate->setValue(values[i].r);
390     vol  ->setValue(values[i].v);
391     option = ext::shared_ptr<VanillaOption>(
392                                         new EuropeanOption(payoff, exercise));
393     option->setPricingEngine(engine);
394     calculated = option->delta();
395     Real error = std::fabs(calculated-values[i].result);
396     Real tolerance = 1e-4;
397     if (error>tolerance)
398         REPORT_FAILURE("delta", payoff, exercise, values[i].s,
399                        values[i].q, values[i].r, today,
400                        values[i].v, values[i].result, calculated,
401                        error, tolerance);
402 
403     i++;
404     payoff = ext::shared_ptr<StrikedTypePayoff>(new
405         PlainVanillaPayoff(values[i].type, values[i].strike));
406     exDate = today + timeToDays(values[i].t);
407     exercise = ext::shared_ptr<Exercise>(new EuropeanExercise(exDate));
408     spot ->setValue(values[i].s);
409     qRate->setValue(values[i].q);
410     rRate->setValue(values[i].r);
411     vol  ->setValue(values[i].v);
412     option = ext::shared_ptr<VanillaOption>(
413                                         new EuropeanOption(payoff, exercise));
414     option->setPricingEngine(engine);
415     calculated = option->delta();
416     error = std::fabs(calculated-values[i].result);
417     if (error>tolerance)
418         REPORT_FAILURE("delta", payoff, exercise, values[i].s,
419                        values[i].q, values[i].r, today,
420                        values[i].v, values[i].result, calculated,
421                        error, tolerance);
422 
423     i++;
424     payoff = ext::shared_ptr<StrikedTypePayoff>(new
425         PlainVanillaPayoff(values[i].type, values[i].strike));
426     exDate = today + timeToDays(values[i].t);
427     exercise = ext::shared_ptr<Exercise>(new EuropeanExercise(exDate));
428     spot ->setValue(values[i].s);
429     qRate->setValue(values[i].q);
430     rRate->setValue(values[i].r);
431     vol  ->setValue(values[i].v);
432     option = ext::shared_ptr<VanillaOption>(
433                                         new EuropeanOption(payoff, exercise));
434     option->setPricingEngine(engine);
435     calculated = option->elasticity();
436     error = std::fabs(calculated-values[i].result);
437     if (error>tolerance)
438         REPORT_FAILURE("elasticity", payoff, exercise, values[i].s,
439                        values[i].q, values[i].r, today,
440                        values[i].v, values[i].result, calculated,
441                        error, tolerance);
442 
443 
444     i++;
445     payoff = ext::shared_ptr<StrikedTypePayoff>(new
446         PlainVanillaPayoff(values[i].type, values[i].strike));
447     exDate = today + timeToDays(values[i].t);
448     exercise = ext::shared_ptr<Exercise>(new EuropeanExercise(exDate));
449     spot ->setValue(values[i].s);
450     qRate->setValue(values[i].q);
451     rRate->setValue(values[i].r);
452     vol  ->setValue(values[i].v);
453     option = ext::shared_ptr<VanillaOption>(
454                                         new EuropeanOption(payoff, exercise));
455     option->setPricingEngine(engine);
456     calculated = option->gamma();
457     error = std::fabs(calculated-values[i].result);
458     if (error>tolerance)
459         REPORT_FAILURE("gamma", payoff, exercise, values[i].s,
460                        values[i].q, values[i].r, today,
461                        values[i].v, values[i].result, calculated,
462                        error, tolerance);
463 
464     i++;
465     payoff = ext::shared_ptr<StrikedTypePayoff>(new
466         PlainVanillaPayoff(values[i].type, values[i].strike));
467     exDate = today + timeToDays(values[i].t);
468     exercise = ext::shared_ptr<Exercise>(new EuropeanExercise(exDate));
469     spot ->setValue(values[i].s);
470     qRate->setValue(values[i].q);
471     rRate->setValue(values[i].r);
472     vol  ->setValue(values[i].v);
473     option = ext::shared_ptr<VanillaOption>(
474                                         new EuropeanOption(payoff, exercise));
475     option->setPricingEngine(engine);
476     calculated = option->gamma();
477     error = std::fabs(calculated-values[i].result);
478     if (error>tolerance)
479         REPORT_FAILURE("gamma", payoff, exercise, values[i].s,
480                        values[i].q, values[i].r, today,
481                        values[i].v, values[i].result, calculated,
482                        error, tolerance);
483 
484 
485     i++;
486     payoff = ext::shared_ptr<StrikedTypePayoff>(new
487         PlainVanillaPayoff(values[i].type, values[i].strike));
488     exDate = today + timeToDays(values[i].t);
489     exercise = ext::shared_ptr<Exercise>(new EuropeanExercise(exDate));
490     spot ->setValue(values[i].s);
491     qRate->setValue(values[i].q);
492     rRate->setValue(values[i].r);
493     vol  ->setValue(values[i].v);
494     option = ext::shared_ptr<VanillaOption>(
495                                         new EuropeanOption(payoff, exercise));
496     option->setPricingEngine(engine);
497     calculated = option->vega();
498     error = std::fabs(calculated-values[i].result);
499     if (error>tolerance)
500         REPORT_FAILURE("vega", payoff, exercise, values[i].s,
501                        values[i].q, values[i].r, today,
502                        values[i].v, values[i].result, calculated,
503                        error, tolerance);
504 
505 
506     i++;
507     payoff = ext::shared_ptr<StrikedTypePayoff>(new
508         PlainVanillaPayoff(values[i].type, values[i].strike));
509     exDate = today + timeToDays(values[i].t);
510     exercise = ext::shared_ptr<Exercise>(new EuropeanExercise(exDate));
511     spot ->setValue(values[i].s);
512     qRate->setValue(values[i].q);
513     rRate->setValue(values[i].r);
514     vol  ->setValue(values[i].v);
515     option = ext::shared_ptr<VanillaOption>(
516                                         new EuropeanOption(payoff, exercise));
517     option->setPricingEngine(engine);
518     calculated = option->vega();
519     error = std::fabs(calculated-values[i].result);
520     if (error>tolerance)
521         REPORT_FAILURE("vega", payoff, exercise, values[i].s,
522                        values[i].q, values[i].r, today,
523                        values[i].v, values[i].result, calculated,
524                        error, tolerance);
525 
526 
527     i++;
528     payoff = ext::shared_ptr<StrikedTypePayoff>(new
529         PlainVanillaPayoff(values[i].type, values[i].strike));
530     exDate = today + timeToDays(values[i].t);
531     exercise = ext::shared_ptr<Exercise>(new EuropeanExercise(exDate));
532     spot ->setValue(values[i].s);
533     qRate->setValue(values[i].q);
534     rRate->setValue(values[i].r);
535     vol  ->setValue(values[i].v);
536     option = ext::shared_ptr<VanillaOption>(
537                                         new EuropeanOption(payoff, exercise));
538     option->setPricingEngine(engine);
539     calculated = option->theta();
540     error = std::fabs(calculated-values[i].result);
541     if (error>tolerance)
542         REPORT_FAILURE("theta", payoff, exercise, values[i].s,
543                        values[i].q, values[i].r, today,
544                        values[i].v, values[i].result, calculated,
545                        error, tolerance);
546 
547 
548     i++;
549     payoff = ext::shared_ptr<StrikedTypePayoff>(new
550         PlainVanillaPayoff(values[i].type, values[i].strike));
551     exDate = today + timeToDays(values[i].t);
552     exercise = ext::shared_ptr<Exercise>(new EuropeanExercise(exDate));
553     spot ->setValue(values[i].s);
554     qRate->setValue(values[i].q);
555     rRate->setValue(values[i].r);
556     vol  ->setValue(values[i].v);
557     option = ext::shared_ptr<VanillaOption>(
558                                         new EuropeanOption(payoff, exercise));
559     option->setPricingEngine(engine);
560     calculated = option->thetaPerDay();
561     error = std::fabs(calculated-values[i].result);
562     if (error>tolerance)
563         REPORT_FAILURE("thetaPerDay", payoff, exercise, values[i].s,
564                        values[i].q, values[i].r, today,
565                        values[i].v, values[i].result, calculated,
566                        error, tolerance);
567 
568 
569     i++;
570     payoff = ext::shared_ptr<StrikedTypePayoff>(new
571         PlainVanillaPayoff(values[i].type, values[i].strike));
572     exDate = today + timeToDays(values[i].t);
573     exercise = ext::shared_ptr<Exercise>(new EuropeanExercise(exDate));
574     spot ->setValue(values[i].s);
575     qRate->setValue(values[i].q);
576     rRate->setValue(values[i].r);
577     vol  ->setValue(values[i].v);
578     option = ext::shared_ptr<VanillaOption>(
579                                         new EuropeanOption(payoff, exercise));
580     option->setPricingEngine(engine);
581     calculated = option->rho();
582     error = std::fabs(calculated-values[i].result);
583     if (error>tolerance)
584         REPORT_FAILURE("rho", payoff, exercise, values[i].s,
585                        values[i].q, values[i].r, today,
586                        values[i].v, values[i].result, calculated,
587                        error, tolerance);
588 
589 
590     i++;
591     payoff = ext::shared_ptr<StrikedTypePayoff>(new
592         PlainVanillaPayoff(values[i].type, values[i].strike));
593     exDate = today + timeToDays(values[i].t);
594     exercise = ext::shared_ptr<Exercise>(new EuropeanExercise(exDate));
595     spot ->setValue(values[i].s);
596     qRate->setValue(values[i].q);
597     rRate->setValue(values[i].r);
598     vol  ->setValue(values[i].v);
599     option = ext::shared_ptr<VanillaOption>(
600                                         new EuropeanOption(payoff, exercise));
601     option->setPricingEngine(engine);
602     calculated = option->dividendRho();
603     error = std::fabs(calculated-values[i].result);
604     if (error>tolerance)
605         REPORT_FAILURE("dividendRho", payoff, exercise, values[i].s,
606                        values[i].q, values[i].r, today,
607                        values[i].v, values[i].result, calculated,
608                        error, tolerance);
609 
610 }
611 
testGreeks()612 void EuropeanOptionTest::testGreeks() {
613 
614     BOOST_TEST_MESSAGE("Testing analytic European option greeks...");
615 
616     using namespace european_option_test;
617 
618     SavedSettings backup;
619 
620     std::map<std::string,Real> calculated, expected, tolerance;
621     tolerance["delta"]  = 1.0e-5;
622     tolerance["gamma"]  = 1.0e-5;
623     tolerance["theta"]  = 1.0e-5;
624     tolerance["rho"]    = 1.0e-5;
625     tolerance["divRho"] = 1.0e-5;
626     tolerance["vega"]   = 1.0e-5;
627 
628     Option::Type types[] = { Option::Call, Option::Put };
629     Real strikes[] = { 50.0, 99.5, 100.0, 100.5, 150.0 };
630     Real underlyings[] = { 100.0 };
631     Rate qRates[] = { 0.04, 0.05, 0.06 };
632     Rate rRates[] = { 0.01, 0.05, 0.15 };
633     Time residualTimes[] = { 1.0, 2.0 };
634     Volatility vols[] = { 0.11, 0.50, 1.20 };
635 
636     DayCounter dc = Actual360();
637     Date today = Date::todaysDate();
638     Settings::instance().evaluationDate() = today;
639 
640     ext::shared_ptr<SimpleQuote> spot(new SimpleQuote(0.0));
641     ext::shared_ptr<SimpleQuote> qRate(new SimpleQuote(0.0));
642     Handle<YieldTermStructure> qTS(flatRate(qRate, dc));
643     ext::shared_ptr<SimpleQuote> rRate(new SimpleQuote(0.0));
644     Handle<YieldTermStructure> rTS(flatRate(rRate, dc));
645     ext::shared_ptr<SimpleQuote> vol(new SimpleQuote(0.0));
646     Handle<BlackVolTermStructure> volTS(flatVol(vol, dc));
647 
648     ext::shared_ptr<StrikedTypePayoff> payoff;
649 
650     for (Size i=0; i<LENGTH(types); i++) {
651       for (Size j=0; j<LENGTH(strikes); j++) {
652         for (Size k=0; k<LENGTH(residualTimes); k++) {
653           Date exDate = today + timeToDays(residualTimes[k]);
654           ext::shared_ptr<Exercise> exercise(new EuropeanExercise(exDate));
655           for (Size kk=0; kk<4; kk++) {
656               // option to check
657               if (kk==0) {
658                   payoff = ext::shared_ptr<StrikedTypePayoff>(new
659                     PlainVanillaPayoff(types[i], strikes[j]));
660               } else if (kk==1) {
661                   payoff = ext::shared_ptr<StrikedTypePayoff>(new
662                     CashOrNothingPayoff(types[i], strikes[j],
663                     100.0));
664               } else if (kk==2) {
665                   payoff = ext::shared_ptr<StrikedTypePayoff>(new
666                     AssetOrNothingPayoff(types[i], strikes[j]));
667               } else if (kk==3) {
668                   payoff = ext::shared_ptr<StrikedTypePayoff>(new
669                     GapPayoff(types[i], strikes[j], 100.0));
670               }
671 
672               ext::shared_ptr<BlackScholesMertonProcess> stochProcess(
673                             new BlackScholesMertonProcess(Handle<Quote>(spot),
674                                                           qTS, rTS, volTS));
675               ext::shared_ptr<PricingEngine> engine(
676                                     new AnalyticEuropeanEngine(stochProcess));
677               EuropeanOption option(payoff, exercise);
678               option.setPricingEngine(engine);
679 
680               for (Size l=0; l<LENGTH(underlyings); l++) {
681                 for (Size m=0; m<LENGTH(qRates); m++) {
682                   for (Size n=0; n<LENGTH(rRates); n++) {
683                     for (Size p=0; p<LENGTH(vols); p++) {
684                       Real u = underlyings[l];
685                       Rate q = qRates[m],
686                            r = rRates[n];
687                       Volatility v = vols[p];
688                       spot->setValue(u);
689                       qRate->setValue(q);
690                       rRate->setValue(r);
691                       vol->setValue(v);
692 
693                       Real value = option.NPV();
694                       calculated["delta"]  = option.delta();
695                       calculated["gamma"]  = option.gamma();
696                       calculated["theta"]  = option.theta();
697                       calculated["rho"]    = option.rho();
698                       calculated["divRho"] = option.dividendRho();
699                       calculated["vega"]   = option.vega();
700 
701                       if (value > spot->value()*1.0e-5) {
702                           // perturb spot and get delta and gamma
703                           Real du = u*1.0e-4;
704                           spot->setValue(u+du);
705                           Real value_p = option.NPV(),
706                                delta_p = option.delta();
707                           spot->setValue(u-du);
708                           Real value_m = option.NPV(),
709                                delta_m = option.delta();
710                           spot->setValue(u);
711                           expected["delta"] = (value_p - value_m)/(2*du);
712                           expected["gamma"] = (delta_p - delta_m)/(2*du);
713 
714                           // perturb rates and get rho and dividend rho
715                           Spread dr = r*1.0e-4;
716                           rRate->setValue(r+dr);
717                           value_p = option.NPV();
718                           rRate->setValue(r-dr);
719                           value_m = option.NPV();
720                           rRate->setValue(r);
721                           expected["rho"] = (value_p - value_m)/(2*dr);
722 
723                           Spread dq = q*1.0e-4;
724                           qRate->setValue(q+dq);
725                           value_p = option.NPV();
726                           qRate->setValue(q-dq);
727                           value_m = option.NPV();
728                           qRate->setValue(q);
729                           expected["divRho"] = (value_p - value_m)/(2*dq);
730 
731                           // perturb volatility and get vega
732                           Volatility dv = v*1.0e-4;
733                           vol->setValue(v+dv);
734                           value_p = option.NPV();
735                           vol->setValue(v-dv);
736                           value_m = option.NPV();
737                           vol->setValue(v);
738                           expected["vega"] = (value_p - value_m)/(2*dv);
739 
740                           // perturb date and get theta
741                           Time dT = dc.yearFraction(today-1, today+1);
742                           Settings::instance().evaluationDate() = today-1;
743                           value_m = option.NPV();
744                           Settings::instance().evaluationDate() = today+1;
745                           value_p = option.NPV();
746                           Settings::instance().evaluationDate() = today;
747                           expected["theta"] = (value_p - value_m)/dT;
748 
749                           // compare
750                           std::map<std::string,Real>::iterator it;
751                           for (it = calculated.begin();
752                                it != calculated.end(); ++it) {
753                               std::string greek = it->first;
754                               Real expct = expected  [greek],
755                                    calcl = calculated[greek],
756                                    tol   = tolerance [greek];
757                               Real error = relativeError(expct,calcl,u);
758                               if (error>tol) {
759                                   REPORT_FAILURE(greek, payoff, exercise,
760                                                  u, q, r, today, v,
761                                                  expct, calcl, error, tol);
762                               }
763                           }
764                       }
765                     }
766                   }
767                 }
768               }
769             }
770         }
771       }
772     }
773 }
774 
testImpliedVol()775 void EuropeanOptionTest::testImpliedVol() {
776 
777     BOOST_TEST_MESSAGE("Testing European option implied volatility...");
778 
779     using namespace european_option_test;
780 
781     SavedSettings backup;
782 
783     Size maxEvaluations = 100;
784     Real tolerance = 1.0e-6;
785 
786     // test options
787     Option::Type types[] = { Option::Call, Option::Put };
788     Real strikes[] = { 90.0, 99.5, 100.0, 100.5, 110.0 };
789     Integer lengths[] = { 36, 180, 360, 1080 };
790 
791     // test data
792     Real underlyings[] = { 90.0, 95.0, 99.9, 100.0, 100.1, 105.0, 110.0 };
793     Rate qRates[] = { 0.01, 0.05, 0.10 };
794     Rate rRates[] = { 0.01, 0.05, 0.10 };
795     Volatility vols[] = { 0.01, 0.20, 0.30, 0.70, 0.90 };
796 
797     DayCounter dc = Actual360();
798     Date today = Date::todaysDate();
799 
800     ext::shared_ptr<SimpleQuote> spot(new SimpleQuote(0.0));
801     ext::shared_ptr<SimpleQuote> vol(new SimpleQuote(0.0));
802     ext::shared_ptr<BlackVolTermStructure> volTS = flatVol(today, vol, dc);
803     ext::shared_ptr<SimpleQuote> qRate(new SimpleQuote(0.0));
804     ext::shared_ptr<YieldTermStructure> qTS = flatRate(today, qRate, dc);
805     ext::shared_ptr<SimpleQuote> rRate(new SimpleQuote(0.0));
806     ext::shared_ptr<YieldTermStructure> rTS = flatRate(today, rRate, dc);
807 
808     for (Size i=0; i<LENGTH(types); i++) {
809       for (Size j=0; j<LENGTH(strikes); j++) {
810         for (Size k=0; k<LENGTH(lengths); k++) {
811           // option to check
812           Date exDate = today + lengths[k];
813           ext::shared_ptr<Exercise> exercise(new EuropeanExercise(exDate));
814           ext::shared_ptr<StrikedTypePayoff> payoff(
815                                 new PlainVanillaPayoff(types[i], strikes[j]));
816           ext::shared_ptr<VanillaOption> option =
817               makeOption(payoff, exercise, spot, qTS, rTS, volTS,
818                          Analytic, Null<Size>(), Null<Size>());
819 
820           ext::shared_ptr<GeneralizedBlackScholesProcess> process =
821               makeProcess(spot, qTS, rTS,volTS);
822 
823           for (Size l=0; l<LENGTH(underlyings); l++) {
824             for (Size m=0; m<LENGTH(qRates); m++) {
825               for (Size n=0; n<LENGTH(rRates); n++) {
826                 for (Size p=0; p<LENGTH(vols); p++) {
827                   Real u = underlyings[l];
828                   Rate q = qRates[m],
829                        r = rRates[n];
830                   Volatility v = vols[p];
831                   spot->setValue(u);
832                   qRate->setValue(q);
833                   rRate->setValue(r);
834                   vol->setValue(v);
835 
836                   Real value = option->NPV();
837                   Volatility implVol = 0.0; // just to remove a warning...
838                   if (value != 0.0) {
839                       // shift guess somehow
840                       vol->setValue(v*0.5);
841                       if (std::fabs(value-option->NPV()) <= 1.0e-12) {
842                           // flat price vs vol --- pointless (and
843                           // numerically unstable) to solve
844                           continue;
845                       }
846                       try {
847                           implVol = option->impliedVolatility(value,
848                                                               process,
849                                                               tolerance,
850                                                               maxEvaluations);
851                       } catch (std::exception& e) {
852                           BOOST_ERROR(
853                               "\nimplied vol calculation failed:" <<
854                               "\n   option:         " << types[i] <<
855                               "\n   strike:         " << strikes[j] <<
856                               "\n   spot value:     " << u <<
857                               "\n   dividend yield: " << io::rate(q) <<
858                               "\n   risk-free rate: " << io::rate(r) <<
859                               "\n   today:          " << today <<
860                               "\n   maturity:       " << exDate <<
861                               "\n   volatility:     " << io::volatility(v) <<
862                               "\n   option value:   " << value <<
863                               "\n" << e.what());
864                       }
865                       if (std::fabs(implVol-v) > tolerance) {
866                           // the difference might not matter
867                           vol->setValue(implVol);
868                           Real value2 = option->NPV();
869                           Real error = relativeError(value,value2,u);
870                           if (error > tolerance) {
871                               BOOST_ERROR(
872                                   types[i] << " option :\n"
873                                   << "    spot value:          " << u << "\n"
874                                   << "    strike:              "
875                                   << strikes[j] << "\n"
876                                   << "    dividend yield:      "
877                                   << io::rate(q) << "\n"
878                                   << "    risk-free rate:      "
879                                   << io::rate(r) << "\n"
880                                   << "    maturity:            "
881                                   << exDate << "\n\n"
882                                   << "    original volatility: "
883                                   << io::volatility(v) << "\n"
884                                   << "    price:               "
885                                   << value << "\n"
886                                   << "    implied volatility:  "
887                                   << io::volatility(implVol)
888                                   << "\n"
889                                   << "    corresponding price: "
890                                   << value2 << "\n"
891                                   << "    error:               " << error);
892                           }
893                       }
894                   }
895                 }
896               }
897             }
898           }
899         }
900       }
901     }
902 }
903 
904 
testImpliedVolContainment()905 void EuropeanOptionTest::testImpliedVolContainment() {
906 
907     BOOST_TEST_MESSAGE("Testing self-containment of "
908                        "implied volatility calculation...");
909 
910     SavedSettings backup;
911 
912     Size maxEvaluations = 100;
913     Real tolerance = 1.0e-6;
914 
915     // test options
916 
917     DayCounter dc = Actual360();
918     Date today = Date::todaysDate();
919 
920     ext::shared_ptr<SimpleQuote> spot(new SimpleQuote(100.0));
921     Handle<Quote> underlying(spot);
922     ext::shared_ptr<SimpleQuote> qRate(new SimpleQuote(0.05));
923     Handle<YieldTermStructure> qTS(flatRate(today, qRate, dc));
924     ext::shared_ptr<SimpleQuote> rRate(new SimpleQuote(0.03));
925     Handle<YieldTermStructure> rTS(flatRate(today, rRate, dc));
926     ext::shared_ptr<SimpleQuote> vol(new SimpleQuote(0.20));
927     Handle<BlackVolTermStructure> volTS(flatVol(today, vol, dc));
928 
929     Date exerciseDate = today + 1*Years;
930     ext::shared_ptr<Exercise> exercise(new EuropeanExercise(exerciseDate));
931     ext::shared_ptr<StrikedTypePayoff> payoff(
932                                  new PlainVanillaPayoff(Option::Call, 100.0));
933 
934     ext::shared_ptr<BlackScholesMertonProcess> process(
935                   new BlackScholesMertonProcess(underlying, qTS, rTS, volTS));
936     ext::shared_ptr<PricingEngine> engine(
937                                         new AnalyticEuropeanEngine(process));
938     // link to the same stochastic process, which shouldn't be changed
939     // by calling methods of either option
940 
941     ext::shared_ptr<VanillaOption> option1(
942                                         new EuropeanOption(payoff, exercise));
943     option1->setPricingEngine(engine);
944     ext::shared_ptr<VanillaOption> option2(
945                                         new EuropeanOption(payoff, exercise));
946     option2->setPricingEngine(engine);
947 
948     // test
949 
950     Real refValue = option2->NPV();
951 
952     Flag f;
953     f.registerWith(option2);
954 
955     option1->impliedVolatility(refValue*1.5, process,
956                                tolerance, maxEvaluations);
957 
958     if (f.isUp())
959         BOOST_ERROR("implied volatility calculation triggered a change "
960                     "in another instrument");
961 
962     option2->recalculate();
963     if (std::fabs(option2->NPV() - refValue) >= 1.0e-8)
964         BOOST_ERROR("implied volatility calculation changed the value "
965                     << "of another instrument: \n"
966                     << std::setprecision(8)
967                     << "previous value: " << refValue << "\n"
968                     << "current value:  " << option2->NPV());
969 
970     vol->setValue(vol->value()*1.5);
971 
972     if (!f.isUp())
973         BOOST_ERROR("volatility change not notified");
974 
975     if (std::fabs(option2->NPV() - refValue) <= 1.0e-8)
976         BOOST_ERROR("volatility change did not cause the value to change");
977 
978 }
979 
980 
981 // different engines
982 
983 namespace {
984 
testEngineConsistency(european_option_test::EngineType engine,Size binomialSteps,Size samples,std::map<std::string,Real> tolerance,bool testGreeks=false)985     void testEngineConsistency(european_option_test::EngineType engine,
986                                Size binomialSteps,
987                                Size samples,
988                                std::map<std::string,Real> tolerance,
989                                bool testGreeks = false) {
990 
991         using namespace european_option_test;
992 
993         std::map<std::string,Real> calculated, expected;
994 
995         // test options
996         Option::Type types[] = { Option::Call, Option::Put };
997         Real strikes[] = { 75.0, 100.0, 125.0 };
998         Integer lengths[] = { 1 };
999 
1000         // test data
1001         Real underlyings[] = { 100.0 };
1002         Rate qRates[] = { 0.00, 0.05 };
1003         Rate rRates[] = { 0.01, 0.05, 0.15 };
1004         Volatility vols[] = { 0.11, 0.50, 1.20 };
1005 
1006         DayCounter dc = Actual360();
1007         Date today = Date::todaysDate();
1008 
1009         ext::shared_ptr<SimpleQuote> spot(new SimpleQuote(0.0));
1010         ext::shared_ptr<SimpleQuote> vol(new SimpleQuote(0.0));
1011         ext::shared_ptr<BlackVolTermStructure> volTS = flatVol(today,vol,dc);
1012         ext::shared_ptr<SimpleQuote> qRate(new SimpleQuote(0.0));
1013         ext::shared_ptr<YieldTermStructure> qTS = flatRate(today,qRate,dc);
1014         ext::shared_ptr<SimpleQuote> rRate(new SimpleQuote(0.0));
1015         ext::shared_ptr<YieldTermStructure> rTS = flatRate(today,rRate,dc);
1016 
1017         for (Size i=0; i<LENGTH(types); i++) {
1018           for (Size j=0; j<LENGTH(strikes); j++) {
1019             for (Size k=0; k<LENGTH(lengths); k++) {
1020               Date exDate = today + lengths[k]*360;
1021               ext::shared_ptr<Exercise> exercise(
1022                                                 new EuropeanExercise(exDate));
1023               ext::shared_ptr<StrikedTypePayoff> payoff(new
1024                                     PlainVanillaPayoff(types[i], strikes[j]));
1025               // reference option
1026               ext::shared_ptr<VanillaOption> refOption =
1027                   makeOption(payoff, exercise, spot, qTS, rTS, volTS,
1028                              Analytic, Null<Size>(), Null<Size>());
1029               // option to check
1030               ext::shared_ptr<VanillaOption> option =
1031                   makeOption(payoff, exercise, spot, qTS, rTS, volTS,
1032                              engine, binomialSteps, samples);
1033 
1034               for (Size l=0; l<LENGTH(underlyings); l++) {
1035                 for (Size m=0; m<LENGTH(qRates); m++) {
1036                   for (Size n=0; n<LENGTH(rRates); n++) {
1037                     for (Size p=0; p<LENGTH(vols); p++) {
1038                       Real u = underlyings[l];
1039                       Rate q = qRates[m],
1040                            r = rRates[n];
1041                       Volatility v = vols[p];
1042                       spot->setValue(u);
1043                       qRate->setValue(q);
1044                       rRate->setValue(r);
1045                       vol->setValue(v);
1046 
1047                       expected.clear();
1048                       calculated.clear();
1049 
1050                       // FLOATING_POINT_EXCEPTION
1051                       expected["value"] = refOption->NPV();
1052                       calculated["value"] = option->NPV();
1053 
1054                       if (testGreeks && option->NPV() > spot->value()*1.0e-5) {
1055                            expected["delta"] = refOption->delta();
1056                            expected["gamma"] = refOption->gamma();
1057                            expected["theta"] = refOption->theta();
1058                            calculated["delta"] = option->delta();
1059                            calculated["gamma"] = option->gamma();
1060                            calculated["theta"] = option->theta();
1061                       }
1062                       std::map<std::string,Real>::iterator it;
1063                       for (it = calculated.begin();
1064                            it != calculated.end(); ++it) {
1065                           std::string greek = it->first;
1066                           Real expct = expected  [greek],
1067                                calcl = calculated[greek],
1068                                tol   = tolerance [greek];
1069                           Real error = relativeError(expct,calcl,u);
1070                           if (error > tol) {
1071                               REPORT_FAILURE(greek, payoff, exercise,
1072                                              u, q, r, today, v,
1073                                              expct, calcl, error, tol);
1074                           }
1075                       }
1076                     }
1077                   }
1078                 }
1079               }
1080             }
1081           }
1082         }
1083     }
1084 
1085 }
1086 
1087 
testJRBinomialEngines()1088 void EuropeanOptionTest::testJRBinomialEngines() {
1089 
1090     BOOST_TEST_MESSAGE("Testing JR binomial European engines "
1091                        "against analytic results...");
1092 
1093     using namespace european_option_test;
1094 
1095     SavedSettings backup;
1096 
1097     EngineType engine = JR;
1098     Size steps = 251;
1099     Size samples = Null<Size>();
1100     std::map<std::string,Real> relativeTol;
1101     relativeTol["value"] = 0.002;
1102     relativeTol["delta"] = 1.0e-3;
1103     relativeTol["gamma"] = 1.0e-4;
1104     relativeTol["theta"] = 0.03;
1105     testEngineConsistency(engine,steps,samples,relativeTol,true);
1106 }
1107 
testCRRBinomialEngines()1108 void EuropeanOptionTest::testCRRBinomialEngines() {
1109 
1110     BOOST_TEST_MESSAGE("Testing CRR binomial European engines "
1111                        "against analytic results...");
1112 
1113     using namespace european_option_test;
1114 
1115     SavedSettings backup;
1116 
1117     EngineType engine = CRR;
1118     Size steps = 501;
1119     Size samples = Null<Size>();
1120     std::map<std::string,Real> relativeTol;
1121     relativeTol["value"] = 0.02;
1122     relativeTol["delta"] = 1.0e-3;
1123     relativeTol["gamma"] = 1.0e-4;
1124     relativeTol["theta"] = 0.03;
1125     testEngineConsistency(engine,steps,samples,relativeTol,true);
1126 }
1127 
testEQPBinomialEngines()1128 void EuropeanOptionTest::testEQPBinomialEngines() {
1129 
1130     BOOST_TEST_MESSAGE("Testing EQP binomial European engines "
1131                        "against analytic results...");
1132 
1133     using namespace european_option_test;
1134 
1135     SavedSettings backup;
1136 
1137     EngineType engine = EQP;
1138     Size steps = 501;
1139     Size samples = Null<Size>();
1140     std::map<std::string,Real> relativeTol;
1141     relativeTol["value"] = 0.02;
1142     relativeTol["delta"] = 1.0e-3;
1143     relativeTol["gamma"] = 1.0e-4;
1144     relativeTol["theta"] = 0.03;
1145     testEngineConsistency(engine,steps,samples,relativeTol,true);
1146 }
1147 
testTGEOBinomialEngines()1148 void EuropeanOptionTest::testTGEOBinomialEngines() {
1149 
1150     BOOST_TEST_MESSAGE("Testing TGEO binomial European engines "
1151                        "against analytic results...");
1152 
1153     using namespace european_option_test;
1154 
1155     SavedSettings backup;
1156 
1157     EngineType engine = TGEO;
1158     Size steps = 251;
1159     Size samples = Null<Size>();
1160     std::map<std::string,Real> relativeTol;
1161     relativeTol["value"] = 0.002;
1162     relativeTol["delta"] = 1.0e-3;
1163     relativeTol["gamma"] = 1.0e-4;
1164     relativeTol["theta"] = 0.03;
1165     testEngineConsistency(engine,steps,samples,relativeTol,true);
1166 }
1167 
testTIANBinomialEngines()1168 void EuropeanOptionTest::testTIANBinomialEngines() {
1169 
1170     BOOST_TEST_MESSAGE("Testing TIAN binomial European engines "
1171                        "against analytic results...");
1172 
1173     using namespace european_option_test;
1174 
1175     SavedSettings backup;
1176 
1177     EngineType engine = TIAN;
1178     Size steps = 251;
1179     Size samples = Null<Size>();
1180     std::map<std::string,Real> relativeTol;
1181     relativeTol["value"] = 0.002;
1182     relativeTol["delta"] = 1.0e-3;
1183     relativeTol["gamma"] = 1.0e-4;
1184     relativeTol["theta"] = 0.03;
1185     testEngineConsistency(engine,steps,samples,relativeTol,true);
1186 }
1187 
testLRBinomialEngines()1188 void EuropeanOptionTest::testLRBinomialEngines() {
1189 
1190     BOOST_TEST_MESSAGE("Testing LR binomial European engines "
1191                        "against analytic results...");
1192 
1193     using namespace european_option_test;
1194 
1195     SavedSettings backup;
1196 
1197     EngineType engine = LR;
1198     Size steps = 251;
1199     Size samples = Null<Size>();
1200     std::map<std::string,Real> relativeTol;
1201     relativeTol["value"] = 1.0e-6;
1202     relativeTol["delta"] = 1.0e-3;
1203     relativeTol["gamma"] = 1.0e-4;
1204     relativeTol["theta"] = 0.03;
1205     testEngineConsistency(engine,steps,samples,relativeTol,true);
1206 }
1207 
testJOSHIBinomialEngines()1208 void EuropeanOptionTest::testJOSHIBinomialEngines() {
1209 
1210     BOOST_TEST_MESSAGE("Testing Joshi binomial European engines "
1211                        "against analytic results...");
1212 
1213     using namespace european_option_test;
1214 
1215     SavedSettings backup;
1216 
1217     EngineType engine = JOSHI;
1218     Size steps = 251;
1219     Size samples = Null<Size>();
1220     std::map<std::string,Real> relativeTol;
1221     relativeTol["value"] = 1.0e-7;
1222     relativeTol["delta"] = 1.0e-3;
1223     relativeTol["gamma"] = 1.0e-4;
1224     relativeTol["theta"] = 0.03;
1225     testEngineConsistency(engine,steps,samples,relativeTol,true);
1226 }
1227 
testFdEngines()1228 void EuropeanOptionTest::testFdEngines() {
1229 
1230     BOOST_TEST_MESSAGE("Testing finite-difference European engines "
1231                        "against analytic results...");
1232 
1233     using namespace european_option_test;
1234 
1235     SavedSettings backup;
1236 
1237     EngineType engine = FiniteDifferences;
1238     Size timeSteps = 500;
1239     Size gridPoints = 500;
1240     std::map<std::string,Real> relativeTol;
1241     relativeTol["value"] = 1.0e-4;
1242     relativeTol["delta"] = 1.0e-6;
1243     relativeTol["gamma"] = 1.0e-6;
1244     relativeTol["theta"] = 1.0e-3;
1245     testEngineConsistency(engine,timeSteps,gridPoints,relativeTol,true);
1246 }
1247 
testIntegralEngines()1248 void EuropeanOptionTest::testIntegralEngines() {
1249 
1250     BOOST_TEST_MESSAGE("Testing integral engines against analytic results...");
1251 
1252     using namespace european_option_test;
1253 
1254     SavedSettings backup;
1255 
1256     EngineType engine = Integral;
1257     Size timeSteps = 300;
1258     Size gridPoints = 300;
1259     std::map<std::string,Real> relativeTol;
1260     relativeTol["value"] = 0.0001;
1261     testEngineConsistency(engine,timeSteps,gridPoints,relativeTol);
1262 }
1263 
testMcEngines()1264 void EuropeanOptionTest::testMcEngines() {
1265 
1266     BOOST_TEST_MESSAGE("Testing Monte Carlo European engines "
1267                        "against analytic results...");
1268 
1269     using namespace european_option_test;
1270 
1271     SavedSettings backup;
1272 
1273     EngineType engine = PseudoMonteCarlo;
1274     Size steps = Null<Size>();
1275     Size samples = 40000;
1276     std::map<std::string,Real> relativeTol;
1277     relativeTol["value"] = 0.01;
1278     testEngineConsistency(engine,steps,samples,relativeTol);
1279 }
1280 
testQmcEngines()1281 void EuropeanOptionTest::testQmcEngines() {
1282 
1283     BOOST_TEST_MESSAGE("Testing Quasi Monte Carlo European engines "
1284                        "against analytic results...");
1285 
1286     using namespace european_option_test;
1287 
1288     SavedSettings backup;
1289 
1290     EngineType engine = QuasiMonteCarlo;
1291     Size steps = Null<Size>();
1292     Size samples = 4095; // 2^12-1
1293     std::map<std::string,Real> relativeTol;
1294     relativeTol["value"] = 0.01;
1295     testEngineConsistency(engine,steps,samples,relativeTol);
1296 }
1297 
testFFTEngines()1298 void EuropeanOptionTest::testFFTEngines() {
1299 
1300     BOOST_TEST_MESSAGE("Testing FFT European engines "
1301                        "against analytic results...");
1302 
1303     using namespace european_option_test;
1304 
1305     SavedSettings backup;
1306 
1307     EngineType engine = FFT;
1308     Size steps = Null<Size>();
1309     Size samples = Null<Size>();
1310     std::map<std::string,Real> relativeTol;
1311     relativeTol["value"] = 0.01;
1312     testEngineConsistency(engine,steps,samples,relativeTol);
1313 }
1314 
1315 
testLocalVolatility()1316 void EuropeanOptionTest::testLocalVolatility() {
1317     BOOST_TEST_MESSAGE("Testing finite-differences with local volatility...");
1318 
1319     using namespace european_option_test;
1320 
1321     SavedSettings backup;
1322 
1323     const Date settlementDate(5, July, 2002);
1324     Settings::instance().evaluationDate() = settlementDate;
1325 
1326     const DayCounter dayCounter = Actual365Fixed();
1327     const Calendar calendar = TARGET();
1328 
1329     Integer t[] = { 13, 41, 75, 165, 256, 345, 524, 703 };
1330     Rate r[] = { 0.0357,0.0349,0.0341,0.0355,0.0359,0.0368,0.0386,0.0401 };
1331 
1332     std::vector<Rate> rates(1, 0.0357);
1333     std::vector<Date> dates(1, settlementDate);
1334     for (Size i = 0; i < 8; ++i) {
1335         dates.push_back(settlementDate + t[i]);
1336         rates.push_back(r[i]);
1337     }
1338     const ext::shared_ptr<YieldTermStructure> rTS(
1339                                    new ZeroCurve(dates, rates, dayCounter));
1340     const ext::shared_ptr<YieldTermStructure> qTS(
1341                                    flatRate(settlementDate, 0.0, dayCounter));
1342 
1343     const ext::shared_ptr<Quote> s0(new SimpleQuote(4500.00));
1344 
1345     Real tmp[] = { 100 ,500 ,2000,3400,3600,3800,4000,4200,4400,4500,
1346                    4600,4800,5000,5200,5400,5600,7500,10000,20000,30000 };
1347     const std::vector<Real> strikes(tmp, tmp+LENGTH(tmp));
1348 
1349     Volatility v[] =
1350       { 1.015873, 1.015873, 1.015873, 0.89729, 0.796493, 0.730914, 0.631335, 0.568895,
1351         0.711309, 0.711309, 0.711309, 0.641309, 0.635593, 0.583653, 0.508045, 0.463182,
1352         0.516034, 0.500534, 0.500534, 0.500534, 0.448706, 0.416661, 0.375470, 0.353442,
1353         0.516034, 0.482263, 0.447713, 0.387703, 0.355064, 0.337438, 0.316966, 0.306859,
1354         0.497587, 0.464373, 0.430764, 0.374052, 0.344336, 0.328607, 0.310619, 0.301865,
1355         0.479511, 0.446815, 0.414194, 0.361010, 0.334204, 0.320301, 0.304664, 0.297180,
1356         0.461866, 0.429645, 0.398092, 0.348638, 0.324680, 0.312512, 0.299082, 0.292785,
1357         0.444801, 0.413014, 0.382634, 0.337026, 0.315788, 0.305239, 0.293855, 0.288660,
1358         0.428604, 0.397219, 0.368109, 0.326282, 0.307555, 0.298483, 0.288972, 0.284791,
1359         0.420971, 0.389782, 0.361317, 0.321274, 0.303697, 0.295302, 0.286655, 0.282948,
1360         0.413749, 0.382754, 0.354917, 0.316532, 0.300016, 0.292251, 0.284420, 0.281164,
1361         0.400889, 0.370272, 0.343525, 0.307904, 0.293204, 0.286549, 0.280189, 0.277767,
1362         0.390685, 0.360399, 0.334344, 0.300507, 0.287149, 0.281380, 0.276271, 0.274588,
1363         0.383477, 0.353434, 0.327580, 0.294408, 0.281867, 0.276746, 0.272655, 0.271617,
1364         0.379106, 0.349214, 0.323160, 0.289618, 0.277362, 0.272641, 0.269332, 0.268846,
1365         0.377073, 0.347258, 0.320776, 0.286077, 0.273617, 0.269057, 0.266293, 0.266265,
1366         0.399925, 0.369232, 0.338895, 0.289042, 0.265509, 0.255589, 0.249308, 0.249665,
1367         0.423432, 0.406891, 0.373720, 0.314667, 0.281009, 0.263281, 0.246451, 0.242166,
1368         0.453704, 0.453704, 0.453704, 0.381255, 0.334578, 0.305527, 0.268909, 0.251367,
1369         0.517748, 0.517748, 0.517748, 0.416577, 0.364770, 0.331595, 0.287423, 0.264285 };
1370 
1371     Matrix blackVolMatrix(strikes.size(), dates.size()-1);
1372     for (Size i=0; i < strikes.size(); ++i)
1373         for (Size j=1; j < dates.size(); ++j) {
1374             blackVolMatrix[i][j-1] = v[i*(dates.size()-1)+j-1];
1375         }
1376 
1377     const ext::shared_ptr<BlackVarianceSurface> volTS(
1378         new BlackVarianceSurface(settlementDate, calendar,
1379                                  std::vector<Date>(dates.begin()+1, dates.end()),
1380                                  strikes, blackVolMatrix,
1381                                  dayCounter));
1382     volTS->setInterpolation<Bicubic>();
1383     const ext::shared_ptr<GeneralizedBlackScholesProcess> process =
1384                                               makeProcess(s0, qTS, rTS,volTS);
1385 
1386     const std::pair<FdmSchemeDesc, std::string> schemeDescs[]= {
1387         std::make_pair(FdmSchemeDesc::Douglas(), "Douglas"),
1388         std::make_pair(FdmSchemeDesc::CrankNicolson(), "Crank-Nicolson"),
1389         std::make_pair(FdmSchemeDesc::ModifiedCraigSneyd(), "Mod. Craig-Sneyd")
1390     };
1391 
1392     for (Size i=2; i < dates.size(); i+=2) {
1393         for (Size j=3; j < strikes.size()-5; j+=5) {
1394             const Date& exDate = dates[i];
1395             const ext::shared_ptr<StrikedTypePayoff> payoff(new
1396                                  PlainVanillaPayoff(Option::Call, strikes[j]));
1397 
1398             const ext::shared_ptr<Exercise> exercise(
1399                                                  new EuropeanExercise(exDate));
1400 
1401             EuropeanOption option(payoff, exercise);
1402             option.setPricingEngine(ext::shared_ptr<PricingEngine>(
1403                                          new AnalyticEuropeanEngine(process)));
1404 
1405             const Real tol = 0.001;
1406             const Real expectedNPV   = option.NPV();
1407             const Real expectedDelta = option.delta();
1408             const Real expectedGamma = option.gamma();
1409 
1410             option.setPricingEngine(ext::shared_ptr<PricingEngine>(
1411                          new FdBlackScholesVanillaEngine(process, 200, 400)));
1412 
1413             Real calculatedNPV = option.NPV();
1414             const Real calculatedDelta = option.delta();
1415             const Real calculatedGamma = option.gamma();
1416 
1417             // check implied pricing first
1418             if (std::fabs(expectedNPV - calculatedNPV) > tol*expectedNPV) {
1419                 BOOST_FAIL("Failed to reproduce option price for "
1420                            << "\n    strike:     " << payoff->strike()
1421                            << "\n    maturity:   " << exDate
1422                            << "\n    calculated: " << calculatedNPV
1423                            << "\n    expected:   " << expectedNPV);
1424             }
1425             if (std::fabs(expectedDelta - calculatedDelta) >tol*expectedDelta) {
1426                 BOOST_FAIL("Failed to reproduce option delta for "
1427                            << "\n    strike:     " << payoff->strike()
1428                            << "\n    maturity:   " << exDate
1429                            << "\n    calculated: " << calculatedDelta
1430                            << "\n    expected:   " << expectedDelta);
1431             }
1432             if (std::fabs(expectedGamma - calculatedGamma) >tol*expectedGamma) {
1433                 BOOST_FAIL("Failed to reproduce option gamma for "
1434                            << "\n    strike:     " << payoff->strike()
1435                            << "\n    maturity:   " << exDate
1436                            << "\n    calculated: " << calculatedGamma
1437                            << "\n    expected:   " << expectedGamma);
1438             }
1439 
1440             // check local vol pricing
1441             // delta/gamma are not the same by definition (model implied greeks)
1442             for (Size i=0; i < LENGTH(schemeDescs); ++i) {
1443                 option.setPricingEngine(
1444                     ext::make_shared<FdBlackScholesVanillaEngine>(
1445                         process, 25, 100, 0, schemeDescs[i].first, true, 0.35));
1446 
1447                 calculatedNPV = option.NPV();
1448                 if (std::fabs(expectedNPV - calculatedNPV) > tol*expectedNPV) {
1449                     BOOST_FAIL("Failed to reproduce local vol option price for "
1450                                << "\n    strike:     " << payoff->strike()
1451                                << "\n    maturity:   " << exDate
1452                                << "\n    calculated: " << calculatedNPV
1453                                << "\n    expected:   " << expectedNPV
1454                                << "\n    scheme:     " << schemeDescs[i].second);
1455                 }
1456             }
1457         }
1458     }
1459 }
1460 
testAnalyticEngineDiscountCurve()1461 void EuropeanOptionTest::testAnalyticEngineDiscountCurve() {
1462     BOOST_TEST_MESSAGE(
1463         "Testing separate discount curve for analytic European engine...");
1464 
1465     SavedSettings backup;
1466 
1467     DayCounter dc = Actual360();
1468     Date today = Date::todaysDate();
1469 
1470     ext::shared_ptr<SimpleQuote> spot(new SimpleQuote(1000.0));
1471     ext::shared_ptr<SimpleQuote> qRate(new SimpleQuote(0.01));
1472     ext::shared_ptr<YieldTermStructure> qTS = flatRate(today, qRate, dc);
1473     ext::shared_ptr<SimpleQuote> rRate(new SimpleQuote(0.015));
1474     ext::shared_ptr<YieldTermStructure> rTS = flatRate(today, rRate, dc);
1475     ext::shared_ptr<SimpleQuote> vol(new SimpleQuote(0.02));
1476     ext::shared_ptr<BlackVolTermStructure> volTS = flatVol(today, vol, dc);
1477     ext::shared_ptr<SimpleQuote> discRate(new SimpleQuote(0.015));
1478     ext::shared_ptr<YieldTermStructure> discTS = flatRate(today, discRate, dc);
1479 
1480     ext::shared_ptr<BlackScholesMertonProcess> stochProcess(new
1481         BlackScholesMertonProcess(Handle<Quote>(spot),
1482             Handle<YieldTermStructure>(qTS),
1483             Handle<YieldTermStructure>(rTS),
1484             Handle<BlackVolTermStructure>(volTS)));
1485     ext::shared_ptr<PricingEngine> engineSingleCurve(
1486         new AnalyticEuropeanEngine(stochProcess));
1487     ext::shared_ptr<PricingEngine> engineMultiCurve(
1488         new AnalyticEuropeanEngine(stochProcess,
1489             Handle<YieldTermStructure>(discTS)));
1490 
1491     ext::shared_ptr<StrikedTypePayoff> payoff(new
1492         PlainVanillaPayoff(Option::Call, 1025.0));
1493     Date exDate = today + Period(1, Years);
1494     ext::shared_ptr<Exercise> exercise(new EuropeanExercise(exDate));
1495     EuropeanOption option(payoff, exercise);
1496     Real npvSingleCurve, npvMultiCurve;
1497     option.setPricingEngine(engineSingleCurve);
1498     npvSingleCurve = option.NPV();
1499     option.setPricingEngine(engineMultiCurve);
1500     npvMultiCurve = option.NPV();
1501     // check that NPV is the same regardless of engine interface
1502     BOOST_CHECK_EQUAL(npvSingleCurve, npvMultiCurve);
1503     // check that NPV changes if discount rate is changed
1504     discRate->setValue(0.023);
1505     npvMultiCurve = option.NPV();
1506     BOOST_CHECK_NE(npvSingleCurve, npvMultiCurve);
1507 }
1508 
1509 
testPDESchemes()1510 void EuropeanOptionTest::testPDESchemes() {
1511     BOOST_TEST_MESSAGE("Testing different PDE schemes "
1512             "to solve Black-Scholes PDEs...");
1513 
1514     SavedSettings backup;
1515 
1516     const DayCounter dc = Actual365Fixed();
1517     const Date today = Date(18, February, 2018);
1518 
1519     Settings::instance().evaluationDate() = today;
1520 
1521     const Handle<Quote> spot(ext::make_shared<SimpleQuote>(100.0));
1522     const Handle<YieldTermStructure> qTS(flatRate(today, 0.06, dc));
1523     const Handle<YieldTermStructure> rTS(flatRate(today, 0.10, dc));
1524     const Handle<BlackVolTermStructure> volTS(flatVol(today, 0.35, dc));
1525 
1526     const Date maturity = today + Period(6, Months);
1527 
1528     const ext::shared_ptr<BlackScholesMertonProcess> process =
1529         ext::make_shared<BlackScholesMertonProcess>(
1530             spot, qTS, rTS, volTS);
1531 
1532     const ext::shared_ptr<PricingEngine> analytic =
1533         ext::make_shared<AnalyticEuropeanEngine>(process);
1534 
1535     // Crank-Nicolson and Douglas scheme are the same in one dimension
1536     const ext::shared_ptr<PricingEngine> douglas =
1537         ext::make_shared<FdBlackScholesVanillaEngine>(
1538             process, 15, 100, 0, FdmSchemeDesc::Douglas());
1539 
1540     const ext::shared_ptr<PricingEngine> crankNicolson =
1541         ext::make_shared<FdBlackScholesVanillaEngine>(
1542             process, 15, 100, 0, FdmSchemeDesc::CrankNicolson());
1543 
1544     const ext::shared_ptr<PricingEngine> implicitEuler =
1545         ext::make_shared<FdBlackScholesVanillaEngine>(
1546             process, 500, 100, 0, FdmSchemeDesc::ImplicitEuler());
1547 
1548     const ext::shared_ptr<PricingEngine> explicitEuler =
1549         ext::make_shared<FdBlackScholesVanillaEngine>(
1550             process, 1000, 100, 0, FdmSchemeDesc::ExplicitEuler());
1551 
1552     const ext::shared_ptr<PricingEngine> methodOfLines =
1553         ext::make_shared<FdBlackScholesVanillaEngine>(
1554             process, 1, 100, 0, FdmSchemeDesc::MethodOfLines());
1555 
1556     const ext::shared_ptr<PricingEngine> hundsdorfer =
1557         ext::make_shared<FdBlackScholesVanillaEngine>(
1558             process, 10, 100, 0, FdmSchemeDesc::Hundsdorfer());
1559 
1560     const ext::shared_ptr<PricingEngine> craigSneyd =
1561         ext::make_shared<FdBlackScholesVanillaEngine>(
1562             process, 10, 100, 0, FdmSchemeDesc::CraigSneyd());
1563 
1564     const ext::shared_ptr<PricingEngine> modCraigSneyd =
1565         ext::make_shared<FdBlackScholesVanillaEngine>(
1566             process, 15, 100, 0, FdmSchemeDesc::ModifiedCraigSneyd());
1567 
1568     const ext::shared_ptr<PricingEngine> trBDF2 =
1569         ext::make_shared<FdBlackScholesVanillaEngine>(
1570             process, 15, 100, 0, FdmSchemeDesc::TrBDF2());
1571 
1572 
1573     const std::pair<ext::shared_ptr<PricingEngine>, std::string> engines[]= {
1574         std::make_pair(douglas, "Douglas"),
1575         std::make_pair(crankNicolson, "Crank-Nicolson"),
1576         std::make_pair(implicitEuler, "Implicit-Euler"),
1577         std::make_pair(explicitEuler, "Explicit-Euler"),
1578         std::make_pair(methodOfLines, "Method-of-Lines"),
1579         std::make_pair(hundsdorfer, "Hundsdorfer"),
1580         std::make_pair(craigSneyd, "Craig-Sneyd"),
1581         std::make_pair(modCraigSneyd, "Modified Craig-Sneyd"),
1582         std::make_pair(trBDF2, "TR-BDF2")
1583     };
1584 
1585     const Size nEngines = LENGTH(engines);
1586 
1587     const ext::shared_ptr<PlainVanillaPayoff> payoff(
1588         ext::make_shared<PlainVanillaPayoff>(Option::Put, spot->value()));
1589 
1590     const ext::shared_ptr<Exercise> exercise(
1591         ext::make_shared<EuropeanExercise>(maturity));
1592 
1593     VanillaOption option(payoff, exercise);
1594 
1595     option.setPricingEngine(analytic);
1596     const Real expected = option.NPV();
1597 
1598     const Real tol = 0.006;
1599     for (Size i=0; i < nEngines; ++i) {
1600         option.setPricingEngine(engines[i].first);
1601         const Real calculated = option.NPV();
1602 
1603         const Real diff = std::fabs(expected - calculated);
1604 
1605         if (diff > tol) {
1606             BOOST_FAIL("Failed to reproduce European option values with the "
1607                     << engines[i].second << " PDE scheme"
1608                        << "\n    calculated: " << calculated
1609                        << "\n    expected:   " << expected
1610                        << "\n    difference: " << diff
1611                        << "\n    tolerance:  " << tol);
1612         }
1613     }
1614 
1615     DividendVanillaOption dividendOption(
1616         payoff, exercise,
1617         std::vector<Date>(1, today + Period(3, Months)),
1618         std::vector<Real>(1, 5.0));
1619 
1620     Array dividendPrices(nEngines);
1621     for (Size i=0; i < nEngines; ++i) {
1622         dividendOption.setPricingEngine(engines[i].first);
1623         dividendPrices[i] = dividendOption.NPV();
1624     }
1625 
1626     const Real expectedDiv = std::accumulate(
1627         dividendPrices.begin(), dividendPrices.end(), 0.0)/nEngines;
1628 
1629     for (Size i=0; i < nEngines; ++i) {
1630         const Real calculated = dividendPrices[i];
1631         const Real diff = std::fabs(expectedDiv - calculated);
1632 
1633         if (diff > tol) {
1634             BOOST_FAIL("Failed to reproduce European option values "
1635                     "with dividend and the "
1636                     << engines[i].second << " PDE scheme"
1637                        << "\n    calculated: " << calculated
1638                        << "\n    expected:   " << expectedDiv
1639                        << "\n    difference: " << diff
1640                        << "\n    tolerance:  " << tol);
1641         }
1642     }
1643 
1644     // make sure that Douglas and Crank-Nicolson are giving the same result
1645     const Size idxDouglas = std::distance(engines,
1646         std::find(engines, engines + LENGTH(engines),
1647             std::make_pair(douglas, std::string("Douglas"))));
1648     const Real douglasNPV = dividendPrices[idxDouglas];
1649 
1650     const Size idxCrankNicolson = std::distance(engines,
1651         std::find(engines, engines + LENGTH(engines),
1652             std::make_pair(crankNicolson, std::string("Crank-Nicolson"))));
1653     const Real crankNicolsonNPV = dividendPrices[idxCrankNicolson];
1654 
1655     const Real schemeTol = 1e-12;
1656     const Real schemeDiff = std::fabs(crankNicolsonNPV - douglasNPV);
1657     if (schemeDiff > schemeTol) {
1658         BOOST_FAIL("Failed to reproduce Douglas scheme option values "
1659                 "with the Crank-Nicolson PDE scheme "
1660                    << "\n    Dougles NPV:        " << douglasNPV
1661                    << "\n    Crank-Nicolson NPV: " << crankNicolsonNPV
1662                    << "\n    difference:         " << schemeDiff
1663                    << "\n    tolerance:          " << schemeTol);
1664     }
1665 }
1666 
testFdEngineWithNonConstantParameters()1667 void EuropeanOptionTest::testFdEngineWithNonConstantParameters() {
1668     BOOST_TEST_MESSAGE("Testing finite-difference European engine "
1669                        "with non-constant parameters...");
1670 
1671     SavedSettings backup;
1672 
1673     Real u = 190.0;
1674     Volatility v = 0.20;
1675 
1676     DayCounter dc = Actual360();
1677     Date today = Settings::instance().evaluationDate();
1678 
1679     ext::shared_ptr<SimpleQuote> spot(new SimpleQuote(u));
1680     ext::shared_ptr<BlackVolTermStructure> volTS = flatVol(today,v,dc);
1681 
1682     std::vector<Date> dates(5);
1683     std::vector<Rate> rates(5);
1684     dates[0] = today;     rates[0] = 0.0;
1685     dates[1] = today+90;  rates[1] = 0.001;
1686     dates[2] = today+180; rates[2] = 0.002;
1687     dates[3] = today+270; rates[3] = 0.005;
1688     dates[4] = today+360; rates[4] = 0.01;
1689     ext::shared_ptr<YieldTermStructure> rTS =
1690         ext::make_shared<ForwardCurve>(dates, rates, dc);
1691     Rate r = rTS->zeroRate(dates[4], dc, Continuous);
1692 
1693     ext::shared_ptr<BlackScholesProcess> process =
1694         ext::make_shared<BlackScholesProcess>(Handle<Quote>(spot),
1695                                               Handle<YieldTermStructure>(rTS),
1696                                               Handle<BlackVolTermStructure>(volTS));
1697 
1698     ext::shared_ptr<Exercise> exercise =
1699         ext::make_shared<EuropeanExercise>(today + 360);
1700     ext::shared_ptr<StrikedTypePayoff> payoff =
1701         ext::make_shared<PlainVanillaPayoff>(Option::Call, 190.0);
1702 
1703     EuropeanOption option(payoff, exercise);
1704 
1705     option.setPricingEngine(ext::make_shared<AnalyticEuropeanEngine>(process));
1706     Real expected = option.NPV();
1707 
1708     Size timeSteps = 200;
1709     Size gridPoints = 201;
1710     option.setPricingEngine(ext::make_shared<FdBlackScholesVanillaEngine>(
1711                               process, timeSteps, gridPoints));
1712     Real calculated = option.NPV();
1713 
1714     Real tolerance = 0.01;
1715     Real error = std::fabs(expected-calculated);
1716     if (error > tolerance) {
1717         REPORT_FAILURE("value", payoff, exercise,
1718                        u, 0.0, r, today, v,
1719                        expected, calculated,
1720                        error, tolerance);
1721     }
1722 }
1723 
testDouglasVsCrankNicolson()1724 void EuropeanOptionTest::testDouglasVsCrankNicolson() {
1725     BOOST_TEST_MESSAGE("Testing Douglas vs Crank-Nicolson scheme "
1726                         "for finite-difference European PDE engines...");
1727 
1728     SavedSettings backup;
1729 
1730     const DayCounter dc = Actual365Fixed();
1731     const Date today = Date(5, October, 2018);
1732 
1733     Settings::instance().evaluationDate() = today;
1734 
1735     const Handle<Quote> spot(ext::make_shared<SimpleQuote>(100.0));
1736     const Handle<YieldTermStructure> qTS(flatRate(today, 0.02, dc));
1737     const Handle<YieldTermStructure> rTS(flatRate(today, 0.075, dc));
1738     const Handle<BlackVolTermStructure> volTS(flatVol(today, 0.25, dc));
1739 
1740     const ext::shared_ptr<BlackScholesMertonProcess> process =
1741         ext::make_shared<BlackScholesMertonProcess>(
1742             spot, qTS, rTS, volTS);
1743 
1744     VanillaOption option(
1745         ext::make_shared<PlainVanillaPayoff>(Option::Put, spot->value()+2),
1746         ext::make_shared<EuropeanExercise>(today + Period(6, Months)));
1747 
1748     option.setPricingEngine(
1749         ext::make_shared<AnalyticEuropeanEngine>(process));
1750 
1751     const Real npv = option.NPV();
1752     const Real schemeTol = 1e-12;
1753     const Real npvTol = 1e-2;
1754 
1755     for (Real theta = 0.2; theta < 0.81; theta+=0.1) {
1756         option.setPricingEngine(
1757             ext::make_shared<FdBlackScholesVanillaEngine>(
1758                 process, 500, 100, 0,
1759                 FdmSchemeDesc(FdmSchemeDesc::CrankNicolsonType, theta, 0.0)));
1760         const Real crankNicolsonNPV = option.NPV();
1761 
1762         const Real npvDiff = std::fabs(crankNicolsonNPV - npv);
1763         if (npvDiff > npvTol) {
1764             BOOST_FAIL("Failed to reproduce european option values "
1765                     "with the Crank-Nicolson PDE scheme "
1766                        << "\n    Analytic NPV:       " << npv
1767                        << "\n    Crank-Nicolson NPV: " << crankNicolsonNPV
1768                        << "\n    theta:              " << theta
1769                        << "\n    difference:         " << npvDiff
1770                        << "\n    tolerance:          " << npvTol);
1771         }
1772 
1773         option.setPricingEngine(
1774             ext::make_shared<FdBlackScholesVanillaEngine>(
1775                 process, 500, 100, 0,
1776                 FdmSchemeDesc(FdmSchemeDesc::DouglasType, theta, 0.0)));
1777         const Real douglasNPV = option.NPV();
1778 
1779         const Real schemeDiff = std::fabs(crankNicolsonNPV - douglasNPV);
1780 
1781         if (schemeDiff > schemeTol) {
1782             BOOST_FAIL("Failed to reproduce Douglas scheme option values "
1783                     "with the Crank-Nicolson PDE scheme "
1784                        << "\n    Dougles NPV:        " << douglasNPV
1785                        << "\n    Crank-Nicolson NPV: " << crankNicolsonNPV
1786                        << "\n    difference:         " << schemeDiff
1787                        << "\n    tolerance:          " << schemeTol);
1788         }
1789     }
1790 }
1791 
suite()1792 test_suite* EuropeanOptionTest::suite() {
1793     test_suite* suite = BOOST_TEST_SUITE("European option tests");
1794     suite->add(QUANTLIB_TEST_CASE(&EuropeanOptionTest::testValues));
1795     suite->add(QUANTLIB_TEST_CASE(&EuropeanOptionTest::testGreekValues));
1796     suite->add(QUANTLIB_TEST_CASE(&EuropeanOptionTest::testGreeks));
1797     suite->add(QUANTLIB_TEST_CASE(&EuropeanOptionTest::testImpliedVol));
1798     suite->add(QUANTLIB_TEST_CASE(
1799                            &EuropeanOptionTest::testImpliedVolContainment));
1800     suite->add(QUANTLIB_TEST_CASE(&EuropeanOptionTest::testJRBinomialEngines));
1801     suite->add(QUANTLIB_TEST_CASE(&EuropeanOptionTest::testCRRBinomialEngines));
1802     suite->add(QUANTLIB_TEST_CASE(&EuropeanOptionTest::testEQPBinomialEngines));
1803     suite->add(QUANTLIB_TEST_CASE(
1804                                &EuropeanOptionTest::testTGEOBinomialEngines));
1805     suite->add(QUANTLIB_TEST_CASE(
1806                                &EuropeanOptionTest::testTIANBinomialEngines));
1807     suite->add(QUANTLIB_TEST_CASE(&EuropeanOptionTest::testLRBinomialEngines));
1808     suite->add(QUANTLIB_TEST_CASE(
1809                               &EuropeanOptionTest::testJOSHIBinomialEngines));
1810     suite->add(QUANTLIB_TEST_CASE(&EuropeanOptionTest::testFdEngines));
1811     suite->add(QUANTLIB_TEST_CASE(&EuropeanOptionTest::testIntegralEngines));
1812     suite->add(QUANTLIB_TEST_CASE(&EuropeanOptionTest::testMcEngines));
1813     suite->add(QUANTLIB_TEST_CASE(&EuropeanOptionTest::testQmcEngines));
1814 
1815     suite->add(QUANTLIB_TEST_CASE(&EuropeanOptionTest::testLocalVolatility));
1816 
1817     suite->add(QUANTLIB_TEST_CASE(
1818                        &EuropeanOptionTest::testAnalyticEngineDiscountCurve));
1819     suite->add(QUANTLIB_TEST_CASE(&EuropeanOptionTest::testPDESchemes));
1820     suite->add(QUANTLIB_TEST_CASE(
1821                  &EuropeanOptionTest::testFdEngineWithNonConstantParameters));
1822     suite->add(QUANTLIB_TEST_CASE(
1823                  &EuropeanOptionTest::testDouglasVsCrankNicolson));
1824 
1825     return suite;
1826 }
1827 
experimental()1828 test_suite* EuropeanOptionTest::experimental() {
1829     test_suite* suite = BOOST_TEST_SUITE("European option experimental tests");
1830     suite->add(QUANTLIB_TEST_CASE(&EuropeanOptionTest::testFFTEngines));
1831     return suite;
1832 }
1833