1	/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3	/*
4	Copyright (C) 2003, 2004 Ferdinando Ametrano
5	Copyright (C) 2005, 2007 StatPro Italia srl
6	Copyright (C) 2005 Joseph Wang
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 "americanoption.hpp"
23	#include "utilities.hpp"
24	#include <ql/any.hpp>
25	#include <ql/time/daycounters/actual360.hpp>
26	#include <ql/time/daycounters/thirty360.hpp>
27	#include <ql/instruments/vanillaoption.hpp>
28	#include <ql/math/functional.hpp>
29	#include <ql/math/randomnumbers/rngtraits.hpp>
30	#include <ql/math/distributions/normaldistribution.hpp>
31	#include <ql/math/integrals/integral.hpp>
32	#include <ql/math/integrals/gausslobattointegral.hpp>
33	#include <ql/math/statistics/incrementalstatistics.hpp>
34	#include <ql/pricingengines/vanilla/analyticeuropeanengine.hpp>
35	#include <ql/pricingengines/vanilla/baroneadesiwhaleyengine.hpp>
36	#include <ql/pricingengines/vanilla/bjerksundstenslandengine.hpp>
37	#include <ql/pricingengines/vanilla/juquadraticengine.hpp>
38	#include <ql/pricingengines/vanilla/fdblackscholesvanillaengine.hpp>
39	#include <ql/pricingengines/vanilla/fdblackscholesshoutengine.hpp>
40	#include <ql/pricingengines/vanilla/qdfpamericanengine.hpp>
41	#include <ql/pricingengines/vanilla/qdplusamericanengine.hpp>
42	#include <ql/termstructures/yield/flatforward.hpp>
43	#include <ql/termstructures/volatility/equityfx/blackconstantvol.hpp>
44	#include <ql/utilities/dataformatters.hpp>
45
46	#include <map>
47
48	using namespace QuantLib;
49	using namespace boost::unit_test_framework;
50
51	#undef REPORT_FAILURE
52	#define REPORT_FAILURE(greekName, payoff, exercise, s, q, r, today, \
53	v, expected, calculated, error, tolerance) \
54	BOOST_ERROR(exerciseTypeToString(exercise) << " " \
55	<< payoff->optionType() << " option with " \
56	<< payoffTypeToString(payoff) << " payoff:\n" \
57	<< " spot value: " << s << "\n" \
58	<< " strike: " << payoff->strike() << "\n" \
59	<< " dividend yield: " << io::rate(q) << "\n" \
60	<< " risk-free rate: " << io::rate(r) << "\n" \
61	<< " reference date: " << today << "\n" \
62	<< " maturity: " << exercise->lastDate() << "\n" \
63	<< " volatility: " << io::volatility(v) << "\n\n" \
64	<< std::fixed << std::setprecision(4) \
65	<< " expected " << greekName << ": " << expected << "\n" \
66	<< " calculated " << greekName << ": " << calculated << "\n"\
67	<< std::scientific \
68	<< " error: " << error << "\n" \
69	<< " tolerance: " << tolerance);
70
71	namespace {
72
73	struct AmericanOptionData {
74	Option::Type type;
75	Real strike;
76	Real s; // spot
77	Rate q; // dividend
78	Rate r; // risk-free rate
79	Time t; // time to maturity
80	Volatility v; // volatility
81	Real result; // expected result
82	};
83
84	}
85
86
87	void AmericanOptionTest::testBaroneAdesiWhaleyValues() {
88
89	BOOST_TEST_MESSAGE("Testing Barone-Adesi and Whaley approximation "
90	"for American options...");
91
92	/* The data below are from
93	"Option pricing formulas", E.G. Haug, McGraw-Hill 1998
94	pag 24
95
96	The following values were replicated only up to the second digit
97	by the VB code provided by Haug, which was used as base for the
98	C++ implementation
99
100	*/
101	AmericanOptionData values[] = {
102	// type, strike, spot, q, r, t, vol, value
103	{ .type: Option::Call, .strike: 100.00, .s: 90.00, .q: 0.10, .r: 0.10, .t: 0.10, .v: 0.15, .result: 0.0206 },
104	{ .type: Option::Call, .strike: 100.00, .s: 100.00, .q: 0.10, .r: 0.10, .t: 0.10, .v: 0.15, .result: 1.8771 },
105	{ .type: Option::Call, .strike: 100.00, .s: 110.00, .q: 0.10, .r: 0.10, .t: 0.10, .v: 0.15, .result: 10.0089 },
106	{ .type: Option::Call, .strike: 100.00, .s: 90.00, .q: 0.10, .r: 0.10, .t: 0.10, .v: 0.25, .result: 0.3159 },
107	{ .type: Option::Call, .strike: 100.00, .s: 100.00, .q: 0.10, .r: 0.10, .t: 0.10, .v: 0.25, .result: 3.1280 },
108	{ .type: Option::Call, .strike: 100.00, .s: 110.00, .q: 0.10, .r: 0.10, .t: 0.10, .v: 0.25, .result: 10.3919 },
109	{ .type: Option::Call, .strike: 100.00, .s: 90.00, .q: 0.10, .r: 0.10, .t: 0.10, .v: 0.35, .result: 0.9495 },
110	{ .type: Option::Call, .strike: 100.00, .s: 100.00, .q: 0.10, .r: 0.10, .t: 0.10, .v: 0.35, .result: 4.3777 },
111	{ .type: Option::Call, .strike: 100.00, .s: 110.00, .q: 0.10, .r: 0.10, .t: 0.10, .v: 0.35, .result: 11.1679 },
112	{ .type: Option::Call, .strike: 100.00, .s: 90.00, .q: 0.10, .r: 0.10, .t: 0.50, .v: 0.15, .result: 0.8208 },
113	{ .type: Option::Call, .strike: 100.00, .s: 100.00, .q: 0.10, .r: 0.10, .t: 0.50, .v: 0.15, .result: 4.0842 },
114	{ .type: Option::Call, .strike: 100.00, .s: 110.00, .q: 0.10, .r: 0.10, .t: 0.50, .v: 0.15, .result: 10.8087 },
115	{ .type: Option::Call, .strike: 100.00, .s: 90.00, .q: 0.10, .r: 0.10, .t: 0.50, .v: 0.25, .result: 2.7437 },
116	{ .type: Option::Call, .strike: 100.00, .s: 100.00, .q: 0.10, .r: 0.10, .t: 0.50, .v: 0.25, .result: 6.8015 },
117	{ .type: Option::Call, .strike: 100.00, .s: 110.00, .q: 0.10, .r: 0.10, .t: 0.50, .v: 0.25, .result: 13.0170 },
118	{ .type: Option::Call, .strike: 100.00, .s: 90.00, .q: 0.10, .r: 0.10, .t: 0.50, .v: 0.35, .result: 5.0063 },
119	{ .type: Option::Call, .strike: 100.00, .s: 100.00, .q: 0.10, .r: 0.10, .t: 0.50, .v: 0.35, .result: 9.5106 },
120	{ .type: Option::Call, .strike: 100.00, .s: 110.00, .q: 0.10, .r: 0.10, .t: 0.50, .v: 0.35, .result: 15.5689 },
121	{ .type: Option::Put, .strike: 100.00, .s: 90.00, .q: 0.10, .r: 0.10, .t: 0.10, .v: 0.15, .result: 10.0000 },
122	{ .type: Option::Put, .strike: 100.00, .s: 100.00, .q: 0.10, .r: 0.10, .t: 0.10, .v: 0.15, .result: 1.8770 },
123	{ .type: Option::Put, .strike: 100.00, .s: 110.00, .q: 0.10, .r: 0.10, .t: 0.10, .v: 0.15, .result: 0.0410 },
124	{ .type: Option::Put, .strike: 100.00, .s: 90.00, .q: 0.10, .r: 0.10, .t: 0.10, .v: 0.25, .result: 10.2533 },
125	{ .type: Option::Put, .strike: 100.00, .s: 100.00, .q: 0.10, .r: 0.10, .t: 0.10, .v: 0.25, .result: 3.1277 },
126	{ .type: Option::Put, .strike: 100.00, .s: 110.00, .q: 0.10, .r: 0.10, .t: 0.10, .v: 0.25, .result: 0.4562 },
127	{ .type: Option::Put, .strike: 100.00, .s: 90.00, .q: 0.10, .r: 0.10, .t: 0.10, .v: 0.35, .result: 10.8787 },
128	{ .type: Option::Put, .strike: 100.00, .s: 100.00, .q: 0.10, .r: 0.10, .t: 0.10, .v: 0.35, .result: 4.3777 },
129	{ .type: Option::Put, .strike: 100.00, .s: 110.00, .q: 0.10, .r: 0.10, .t: 0.10, .v: 0.35, .result: 1.2402 },
130	{ .type: Option::Put, .strike: 100.00, .s: 90.00, .q: 0.10, .r: 0.10, .t: 0.50, .v: 0.15, .result: 10.5595 },
131	{ .type: Option::Put, .strike: 100.00, .s: 100.00, .q: 0.10, .r: 0.10, .t: 0.50, .v: 0.15, .result: 4.0842 },
132	{ .type: Option::Put, .strike: 100.00, .s: 110.00, .q: 0.10, .r: 0.10, .t: 0.50, .v: 0.15, .result: 1.0822 },
133	{ .type: Option::Put, .strike: 100.00, .s: 90.00, .q: 0.10, .r: 0.10, .t: 0.50, .v: 0.25, .result: 12.4419 },
134	{ .type: Option::Put, .strike: 100.00, .s: 100.00, .q: 0.10, .r: 0.10, .t: 0.50, .v: 0.25, .result: 6.8014 },
135	{ .type: Option::Put, .strike: 100.00, .s: 110.00, .q: 0.10, .r: 0.10, .t: 0.50, .v: 0.25, .result: 3.3226 },
136	{ .type: Option::Put, .strike: 100.00, .s: 90.00, .q: 0.10, .r: 0.10, .t: 0.50, .v: 0.35, .result: 14.6945 },
137	{ .type: Option::Put, .strike: 100.00, .s: 100.00, .q: 0.10, .r: 0.10, .t: 0.50, .v: 0.35, .result: 9.5104 },
138	{ .type: Option::Put, .strike: 100.00, .s: 110.00, .q: 0.10, .r: 0.10, .t: 0.50, .v: 0.35, .result: 5.8823 },
139	{ .type: Option::Put, .strike: 100.00, .s: 100.00, .q: 0.00, .r: 0.00, .t: 0.50, .v: 0.15, .result: 4.2294 }
140	};
141
142	Date today = Date::todaysDate();
143	DayCounter dc = Actual360();
144	ext::shared_ptr<SimpleQuote> spot(new SimpleQuote(0.0));
145	ext::shared_ptr<SimpleQuote> qRate(new SimpleQuote(0.0));
146	ext::shared_ptr<YieldTermStructure> qTS = flatRate(today, forward: qRate, dc);
147	ext::shared_ptr<SimpleQuote> rRate(new SimpleQuote(0.0));
148	ext::shared_ptr<YieldTermStructure> rTS = flatRate(today, forward: rRate, dc);
149	ext::shared_ptr<SimpleQuote> vol(new SimpleQuote(0.0));
150	ext::shared_ptr<BlackVolTermStructure> volTS = flatVol(today, volatility: vol, dc);
151
152	Real tolerance = 3.0e-3;
153
154	for (auto& value : values) {
155
156	ext::shared_ptr<StrikedTypePayoff> payoff(new PlainVanillaPayoff(value.type, value.strike));
157	Date exDate = today + timeToDays(t: value.t);
158	ext::shared_ptr<Exercise> exercise(
159	new AmericanExercise(today, exDate));
160
161	spot->setValue(value.s);
162	qRate->setValue(value.q);
163	rRate->setValue(value.r);
164	vol->setValue(value.v);
165
166	ext::shared_ptr<BlackScholesMertonProcess> stochProcess(new
167	BlackScholesMertonProcess(Handle<Quote>(spot),
168	Handle<YieldTermStructure>(qTS),
169	Handle<YieldTermStructure>(rTS),
170	Handle<BlackVolTermStructure>(volTS)));
171
172	ext::shared_ptr<PricingEngine> engine(
173	new BaroneAdesiWhaleyApproximationEngine(stochProcess));
174
175	VanillaOption option(payoff, exercise);
176	option.setPricingEngine(engine);
177
178	Real calculated = option.NPV();
179	Real error = std::fabs(x: calculated - value.result);
180	if (error > tolerance) {
181	REPORT_FAILURE("value", payoff, exercise, value.s, value.q, value.r, today, value.v,
182	value.result, calculated, error, tolerance);
183	}
184	}
185	}
186
187
188	void AmericanOptionTest::testBjerksundStenslandValues() {
189
190	BOOST_TEST_MESSAGE("Testing Bjerksund and Stensland approximation "
191	"for American options...");
192
193	AmericanOptionData values[] = {
194	// type, strike, spot, q, r, t, vol, value, tol
195	// from "Option pricing formulas", Haug, McGraw-Hill 1998, pag 27
196	{ .type: Option::Call, .strike: 40.00, .s: 42.00, .q: 0.08, .r: 0.04, .t: 0.75, .v: 0.35, .result: 5.2704 },
197	// from "Option pricing formulas", Haug, McGraw-Hill 1998, VBA code
198	{ .type: Option::Put, .strike: 40.00, .s: 36.00, .q: 0.00, .r: 0.06, .t: 1.00, .v: 0.20, .result: 4.4531 },
199	// ATM option with very small volatility, reference value taken from R
200	{ .type: Option::Call, .strike: 100, .s: 100, .q: 0.05, .r: 0.05, .t: 1.0, .v: 0.0021, .result: 0.08032314 },
201	// ATM option with very small volatility,
202	// reference value taken from Barone-Adesi and Whaley Approximation
203	{ .type: Option::Call, .strike: 100, .s: 100, .q: 0.05, .r: 0.05, .t: 1.0, .v: 0.0001, .result: 0.003860656 },
204	{ .type: Option::Call, .strike: 100, .s: 99.99, .q: 0.05, .r: 0.05, .t: 1.0, .v: 0.0001, .result: 0.00081 },
205	// ITM option with a very small volatility
206	{ .type: Option::Call, .strike: 100, .s: 110, .q: 0.05, .r: 0.05, .t: 1.0, .v: 0.0001, .result: 10.0 },
207	{ .type: Option::Put, .strike: 110, .s: 100, .q: 0.05, .r: 0.05, .t: 1.0, .v: 0.0001, .result: 10.0 },
208	// ATM option with a very large volatility
209	{ .type: Option::Put, .strike: 100, .s: 110, .q: 0.05, .r: 0.05, .t: 1.0, .v: 10, .result: 95.12289 }
210	};
211
212	Date today = Date::todaysDate();
213	DayCounter dc = Actual360();
214	ext::shared_ptr<SimpleQuote> spot(new SimpleQuote(0.0));
215	ext::shared_ptr<SimpleQuote> qRate(new SimpleQuote(0.0));
216	ext::shared_ptr<YieldTermStructure> qTS = flatRate(today, forward: qRate, dc);
217	ext::shared_ptr<SimpleQuote> rRate(new SimpleQuote(0.0));
218	ext::shared_ptr<YieldTermStructure> rTS = flatRate(today, forward: rRate, dc);
219	ext::shared_ptr<SimpleQuote> vol(new SimpleQuote(0.0));
220	ext::shared_ptr<BlackVolTermStructure> volTS = flatVol(today, volatility: vol, dc);
221
222	Real tolerance = 5.0e-5;
223
224	for (auto& value : values) {
225
226	ext::shared_ptr<StrikedTypePayoff> payoff(new PlainVanillaPayoff(value.type, value.strike));
227	Date exDate = today + timeToDays(t: value.t);
228	ext::shared_ptr<Exercise> exercise(new AmericanExercise(today, exDate));
229
230	spot->setValue(value.s);
231	qRate->setValue(value.q);
232	rRate->setValue(value.r);
233	vol->setValue(value.v);
234
235	ext::shared_ptr<BlackScholesMertonProcess> stochProcess(new
236	BlackScholesMertonProcess(Handle<Quote>(spot),
237	Handle<YieldTermStructure>(qTS),
238	Handle<YieldTermStructure>(rTS),
239	Handle<BlackVolTermStructure>(volTS)));
240
241	ext::shared_ptr<PricingEngine> engine(
242	new BjerksundStenslandApproximationEngine(stochProcess));
243
244	VanillaOption option(payoff, exercise);
245	option.setPricingEngine(engine);
246
247	Real calculated = option.NPV();
248	Real error = std::fabs(x: calculated - value.result);
249	if (error > tolerance) {
250	REPORT_FAILURE("value", payoff, exercise, value.s, value.q, value.r, today, value.v,
251	value.result, calculated, error, tolerance);
252	}
253	}
254	}
255
256	namespace {
257
258	/* The data below are from
259	An Approximate Formula for Pricing American Options
260	Journal of Derivatives Winter 1999
261	Ju, N.
262	*/
263	AmericanOptionData juValues[] = {
264	// type, strike, spot, q, r, t, vol, value, tol
265	// These values are from Exhibit 3 - Short dated Put Options
266	{ .type: Option::Put, .strike: 35.00, .s: 40.00, .q: 0.0, .r: 0.0488, .t: 0.0833, .v: 0.2, .result: 0.006 },
267	{ .type: Option::Put, .strike: 35.00, .s: 40.00, .q: 0.0, .r: 0.0488, .t: 0.3333, .v: 0.2, .result: 0.201 },
268	{ .type: Option::Put, .strike: 35.00, .s: 40.00, .q: 0.0, .r: 0.0488, .t: 0.5833, .v: 0.2, .result: 0.433 },
269
270	{ .type: Option::Put, .strike: 40.00, .s: 40.00, .q: 0.0, .r: 0.0488, .t: 0.0833, .v: 0.2, .result: 0.851 },
271	{ .type: Option::Put, .strike: 40.00, .s: 40.00, .q: 0.0, .r: 0.0488, .t: 0.3333, .v: 0.2, .result: 1.576 },
272	{ .type: Option::Put, .strike: 40.00, .s: 40.00, .q: 0.0, .r: 0.0488, .t: 0.5833, .v: 0.2, .result: 1.984 },
273
274	{ .type: Option::Put, .strike: 45.00, .s: 40.00, .q: 0.0, .r: 0.0488, .t: 0.0833, .v: 0.2, .result: 5.000 },
275	{ .type: Option::Put, .strike: 45.00, .s: 40.00, .q: 0.0, .r: 0.0488, .t: 0.3333, .v: 0.2, .result: 5.084 },
276	{ .type: Option::Put, .strike: 45.00, .s: 40.00, .q: 0.0, .r: 0.0488, .t: 0.5833, .v: 0.2, .result: 5.260 },
277
278	{ .type: Option::Put, .strike: 35.00, .s: 40.00, .q: 0.0, .r: 0.0488, .t: 0.0833, .v: 0.3, .result: 0.078 },
279	{ .type: Option::Put, .strike: 35.00, .s: 40.00, .q: 0.0, .r: 0.0488, .t: 0.3333, .v: 0.3, .result: 0.697 },
280	{ .type: Option::Put, .strike: 35.00, .s: 40.00, .q: 0.0, .r: 0.0488, .t: 0.5833, .v: 0.3, .result: 1.218 },
281
282	{ .type: Option::Put, .strike: 40.00, .s: 40.00, .q: 0.0, .r: 0.0488, .t: 0.0833, .v: 0.3, .result: 1.309 },
283	{ .type: Option::Put, .strike: 40.00, .s: 40.00, .q: 0.0, .r: 0.0488, .t: 0.3333, .v: 0.3, .result: 2.477 },
284	{ .type: Option::Put, .strike: 40.00, .s: 40.00, .q: 0.0, .r: 0.0488, .t: 0.5833, .v: 0.3, .result: 3.161 },
285
286	{ .type: Option::Put, .strike: 45.00, .s: 40.00, .q: 0.0, .r: 0.0488, .t: 0.0833, .v: 0.3, .result: 5.059 },
287	{ .type: Option::Put, .strike: 45.00, .s: 40.00, .q: 0.0, .r: 0.0488, .t: 0.3333, .v: 0.3, .result: 5.699 },
288	{ .type: Option::Put, .strike: 45.00, .s: 40.00, .q: 0.0, .r: 0.0488, .t: 0.5833, .v: 0.3, .result: 6.231 },
289
290	{ .type: Option::Put, .strike: 35.00, .s: 40.00, .q: 0.0, .r: 0.0488, .t: 0.0833, .v: 0.4, .result: 0.247 },
291	{ .type: Option::Put, .strike: 35.00, .s: 40.00, .q: 0.0, .r: 0.0488, .t: 0.3333, .v: 0.4, .result: 1.344 },
292	{ .type: Option::Put, .strike: 35.00, .s: 40.00, .q: 0.0, .r: 0.0488, .t: 0.5833, .v: 0.4, .result: 2.150 },
293
294	{ .type: Option::Put, .strike: 40.00, .s: 40.00, .q: 0.0, .r: 0.0488, .t: 0.0833, .v: 0.4, .result: 1.767 },
295	{ .type: Option::Put, .strike: 40.00, .s: 40.00, .q: 0.0, .r: 0.0488, .t: 0.3333, .v: 0.4, .result: 3.381 },
296	{ .type: Option::Put, .strike: 40.00, .s: 40.00, .q: 0.0, .r: 0.0488, .t: 0.5833, .v: 0.4, .result: 4.342 },
297
298	{ .type: Option::Put, .strike: 45.00, .s: 40.00, .q: 0.0, .r: 0.0488, .t: 0.0833, .v: 0.4, .result: 5.288 },
299	{ .type: Option::Put, .strike: 45.00, .s: 40.00, .q: 0.0, .r: 0.0488, .t: 0.3333, .v: 0.4, .result: 6.501 },
300	{ .type: Option::Put, .strike: 45.00, .s: 40.00, .q: 0.0, .r: 0.0488, .t: 0.5833, .v: 0.4, .result: 7.367 },
301
302	// Type in Exhibits 4 and 5 if you have some spare time ;-)
303
304	// type, strike, spot, q, r, t, vol, value, tol
305	// values from Exhibit 6 - Long dated Call Options with dividends
306	{ .type: Option::Call, .strike: 100.00, .s: 80.00, .q: 0.07, .r: 0.03, .t: 3.0, .v: 0.2, .result: 2.605 },
307	{ .type: Option::Call, .strike: 100.00, .s: 90.00, .q: 0.07, .r: 0.03, .t: 3.0, .v: 0.2, .result: 5.182 },
308	{ .type: Option::Call, .strike: 100.00, .s: 100.00, .q: 0.07, .r: 0.03, .t: 3.0, .v: 0.2, .result: 9.065 },
309	{ .type: Option::Call, .strike: 100.00, .s: 110.00, .q: 0.07, .r: 0.03, .t: 3.0, .v: 0.2, .result: 14.430 },
310	{ .type: Option::Call, .strike: 100.00, .s: 120.00, .q: 0.07, .r: 0.03, .t: 3.0, .v: 0.2, .result: 21.398 },
311
312	{ .type: Option::Call, .strike: 100.00, .s: 80.00, .q: 0.07, .r: 0.03, .t: 3.0, .v: 0.4, .result: 11.336 },
313	{ .type: Option::Call, .strike: 100.00, .s: 90.00, .q: 0.07, .r: 0.03, .t: 3.0, .v: 0.4, .result: 15.711 },
314	{ .type: Option::Call, .strike: 100.00, .s: 100.00, .q: 0.07, .r: 0.03, .t: 3.0, .v: 0.4, .result: 20.760 },
315	{ .type: Option::Call, .strike: 100.00, .s: 110.00, .q: 0.07, .r: 0.03, .t: 3.0, .v: 0.4, .result: 26.440 },
316	{ .type: Option::Call, .strike: 100.00, .s: 120.00, .q: 0.07, .r: 0.03, .t: 3.0, .v: 0.4, .result: 32.709 },
317
318	{ .type: Option::Call, .strike: 100.00, .s: 80.00, .q: 0.07, .r: 0.00001, .t: 3.0, .v: 0.3, .result: 5.552 },
319	{ .type: Option::Call, .strike: 100.00, .s: 90.00, .q: 0.07, .r: 0.00001, .t: 3.0, .v: 0.3, .result: 8.868 },
320	{ .type: Option::Call, .strike: 100.00, .s: 100.00, .q: 0.07, .r: 0.00001, .t: 3.0, .v: 0.3, .result: 13.158 },
321	{ .type: Option::Call, .strike: 100.00, .s: 110.00, .q: 0.07, .r: 0.00001, .t: 3.0, .v: 0.3, .result: 18.458 },
322	{ .type: Option::Call, .strike: 100.00, .s: 120.00, .q: 0.07, .r: 0.00001, .t: 3.0, .v: 0.3, .result: 24.786 },
323
324	{ .type: Option::Call, .strike: 100.00, .s: 80.00, .q: 0.03, .r: 0.07, .t: 3.0, .v: 0.3, .result: 12.177 },
325	{ .type: Option::Call, .strike: 100.00, .s: 90.00, .q: 0.03, .r: 0.07, .t: 3.0, .v: 0.3, .result: 17.411 },
326	{ .type: Option::Call, .strike: 100.00, .s: 100.00, .q: 0.03, .r: 0.07, .t: 3.0, .v: 0.3, .result: 23.402 },
327	{ .type: Option::Call, .strike: 100.00, .s: 110.00, .q: 0.03, .r: 0.07, .t: 3.0, .v: 0.3, .result: 30.028 },
328	{ .type: Option::Call, .strike: 100.00, .s: 120.00, .q: 0.03, .r: 0.07, .t: 3.0, .v: 0.3, .result: 37.177 }
329	};
330
331	}
332
333
334	void AmericanOptionTest::testJuValues() {
335
336	BOOST_TEST_MESSAGE("Testing Ju approximation for American options...");
337
338	Date today = Date::todaysDate();
339	DayCounter dc = Actual360();
340	ext::shared_ptr<SimpleQuote> spot(new SimpleQuote(0.0));
341	ext::shared_ptr<SimpleQuote> qRate(new SimpleQuote(0.0));
342	ext::shared_ptr<YieldTermStructure> qTS = flatRate(today, forward: qRate, dc);
343	ext::shared_ptr<SimpleQuote> rRate(new SimpleQuote(0.0));
344	ext::shared_ptr<YieldTermStructure> rTS = flatRate(today, forward: rRate, dc);
345	ext::shared_ptr<SimpleQuote> vol(new SimpleQuote(0.0));
346	ext::shared_ptr<BlackVolTermStructure> volTS = flatVol(today, volatility: vol, dc);
347
348	Real tolerance = 1.0e-3;
349
350	for (auto& juValue : juValues) {
351
352	ext::shared_ptr<StrikedTypePayoff> payoff(
353	new PlainVanillaPayoff(juValue.type, juValue.strike));
354	Date exDate = today + timeToDays(t: juValue.t);
355	ext::shared_ptr<Exercise> exercise(
356	new AmericanExercise(today, exDate));
357
358	spot->setValue(juValue.s);
359	qRate->setValue(juValue.q);
360	rRate->setValue(juValue.r);
361	vol->setValue(juValue.v);
362
363	ext::shared_ptr<BlackScholesMertonProcess> stochProcess(new
364	BlackScholesMertonProcess(Handle<Quote>(spot),
365	Handle<YieldTermStructure>(qTS),
366	Handle<YieldTermStructure>(rTS),
367	Handle<BlackVolTermStructure>(volTS)));
368
369	ext::shared_ptr<PricingEngine> engine(
370	new JuQuadraticApproximationEngine(stochProcess));
371
372	VanillaOption option(payoff, exercise);
373	option.setPricingEngine(engine);
374
375	Real calculated = option.NPV();
376	Real error = std::fabs(x: calculated - juValue.result);
377	if (error > tolerance) {
378	REPORT_FAILURE("value", payoff, exercise, juValue.s, juValue.q, juValue.r, today,
379	juValue.v, juValue.result, calculated, error, tolerance);
380	}
381	}
382	}
383
384
385	void AmericanOptionTest::testFdValues() {
386
387	BOOST_TEST_MESSAGE("Testing finite-difference and QR+ engine "
388	"for American options...");
389
390	Date today = Date::todaysDate();
391	DayCounter dc = Actual360();
392	ext::shared_ptr<SimpleQuote> spot(new SimpleQuote(0.0));
393	ext::shared_ptr<SimpleQuote> qRate(new SimpleQuote(0.0));
394	ext::shared_ptr<YieldTermStructure> qTS = flatRate(today, forward: qRate, dc);
395	ext::shared_ptr<SimpleQuote> rRate(new SimpleQuote(0.0));
396	ext::shared_ptr<YieldTermStructure> rTS = flatRate(today, forward: rRate, dc);
397	ext::shared_ptr<SimpleQuote> vol(new SimpleQuote(0.0));
398	ext::shared_ptr<BlackVolTermStructure> volTS = flatVol(today, volatility: vol, dc);
399
400	Real tolerance = 8.0e-2;
401
402	ext::shared_ptr<BlackScholesMertonProcess> stochProcess =
403	ext::make_shared<BlackScholesMertonProcess>(
404	args: Handle<Quote>(spot),
405	args: Handle<YieldTermStructure>(qTS),
406	args: Handle<YieldTermStructure>(rTS),
407	args: Handle<BlackVolTermStructure>(volTS));
408
409	ext::shared_ptr<PricingEngine> pdeEngine =
410	ext::make_shared<FdBlackScholesVanillaEngine>(args&: stochProcess, args: 100, args: 400);
411
412	ext::shared_ptr<PricingEngine> qrPlusEngine =
413	ext::make_shared<FdBlackScholesVanillaEngine>(args&: stochProcess);
414
415	for (auto& juValue : juValues) {
416
417	ext::shared_ptr<StrikedTypePayoff> payoff(
418	new PlainVanillaPayoff(juValue.type, juValue.strike));
419
420	Date exDate = today + timeToDays(t: juValue.t);
421	ext::shared_ptr<Exercise> exercise(
422	new AmericanExercise(today, exDate));
423
424	spot->setValue(juValue.s);
425	qRate->setValue(juValue.q);
426	rRate->setValue(juValue.r);
427	vol->setValue(juValue.v);
428
429	VanillaOption option(payoff, exercise);
430	option.setPricingEngine(pdeEngine);
431
432	Real pdeCalculated = option.NPV();
433	Real error = std::fabs(x: pdeCalculated - juValue.result);
434	if (error > tolerance) {
435	REPORT_FAILURE("value", payoff, exercise, juValue.s, juValue.q, juValue.r, today,
436	juValue.v, juValue.result, pdeCalculated, error, tolerance);
437	}
438
439	option.setPricingEngine(qrPlusEngine);
440
441	Real qrPlusCalculated = option.NPV();
442	if (std::abs(x: pdeCalculated - qrPlusCalculated) > 2e-2)
443	BOOST_FAIL("QR+ boundary approximation failed to "
444	"reproduce PDE value for "
445	<< "\n OptionType: " <<
446	((juValue.type == Option::Call)? "Call" : "Put")
447	<< std::setprecision(16)
448	<< "\n spot: " << spot->value()
449	<< "\n strike: " << juValue.strike
450	<< "\n r: " << rRate->value()
451	<< "\n q: " << qRate->value()
452	<< "\n vol: " << vol->value()
453	<< "\n PDE value: " << pdeCalculated
454	<< "\n QR+ value: " << qrPlusCalculated);
455	}
456	}
457
458
459	namespace {
460
461	template <class Engine>
462	void testFdGreeks() {
463
464	std::map<std::string,Real> calculated, expected, tolerance;
465	tolerance["delta"] = 7.0e-4;
466	tolerance["gamma"] = 2.0e-4;
467	//tolerance["theta"] = 1.0e-4;
468
469	Option::Type types[] = { Option::Call, Option::Put };
470	Real strikes[] = { 50.0, 99.5, 100.0, 100.5, 150.0 };
471	Real underlyings[] = { 100.0 };
472	Rate qRates[] = { 0.04, 0.05, 0.06 };
473	Rate rRates[] = { 0.01, 0.05, 0.15 };
474	Integer years[] = { 1, 2 };
475	Volatility vols[] = { 0.11, 0.50, 1.20 };
476
477	DayCounter dc = Actual360();
478	Date today = Date::todaysDate();
479	Settings::instance().evaluationDate() = today;
480
481	ext::shared_ptr<SimpleQuote> spot(new SimpleQuote(0.0));
482	ext::shared_ptr<SimpleQuote> qRate(new SimpleQuote(0.0));
483	Handle<YieldTermStructure> qTS(flatRate(forward: qRate, dc));
484	ext::shared_ptr<SimpleQuote> rRate(new SimpleQuote(0.0));
485	Handle<YieldTermStructure> rTS(flatRate(forward: rRate, dc));
486	ext::shared_ptr<SimpleQuote> vol(new SimpleQuote(0.0));
487	Handle<BlackVolTermStructure> volTS(flatVol(volatility: vol, dc));
488
489	ext::shared_ptr<StrikedTypePayoff> payoff;
490
491	for (auto& type : types) {
492	for (Real strike : strikes) {
493	for (int year : years) {
494	Date exDate = today + year * Years;
495	ext::shared_ptr<Exercise> exercise(new AmericanExercise(today, exDate));
496	ext::shared_ptr<StrikedTypePayoff> payoff(new PlainVanillaPayoff(type, strike));
497	ext::shared_ptr<BlackScholesMertonProcess> stochProcess(
498	new BlackScholesMertonProcess(Handle<Quote>(spot), qTS, rTS, volTS));
499
500	ext::shared_ptr<PricingEngine> engine(new Engine(stochProcess, 50));
501
502	VanillaOption option(payoff, exercise);
503	option.setPricingEngine(engine);
504
505	for (Real u : underlyings) {
506	for (Real m : qRates) {
507	for (Real n : rRates) {
508	for (Real v : vols) {
509	Rate q = m, r = n;
510	spot->setValue(u);
511	qRate->setValue(q);
512	rRate->setValue(r);
513	vol->setValue(v);
514	Real value = option.NPV();
515	calculated["delta"] = option.delta();
516	calculated["gamma"] = option.gamma();
517	// calculated["theta"] = option.theta();
518
519	if (value > spot->value() * 1.0e-5) {
520	// perturb spot and get delta and gamma
521	Real du = u * 1.0e-4;
522	spot->setValue(u + du);
523	Real value_p = option.NPV(), delta_p = option.delta();
524	spot->setValue(u - du);
525	Real value_m = option.NPV(), delta_m = option.delta();
526	spot->setValue(u);
527	expected["delta"] = (value_p - value_m) / (2 * du);
528	expected["gamma"] = (delta_p - delta_m) / (2 * du);
529
530	/*
531	// perturb date and get theta
532	Time dT = dc.yearFraction(today-1, today+1);
533	Settings::instance().setEvaluationDate(today-1);
534	value_m = option.NPV();
535	Settings::instance().setEvaluationDate(today+1);
536	value_p = option.NPV();
537	Settings::instance().setEvaluationDate(today);
538	expected["theta"] = (value_p - value_m)/dT;
539	*/
540
541	// compare
542	std::map<std::string, Real>::iterator it;
543	for (it = calculated.begin(); it != calculated.end();
544	++it) {
545	std::string greek = it->first;
546	Real expct = expected[greek], calcl = calculated[greek],
547	tol = tolerance[greek];
548	Real error = relativeError(x1: expct, x2: calcl, reference: u);
549	if (error > tol) {
550	REPORT_FAILURE(greek, payoff, exercise, u, q, r,
551	today, v, expct, calcl, error, tol);
552	}
553	}
554	}
555	}
556	}
557	}
558	}
559	}
560	}
561	}
562	}
563
564	}
565
566
567	void AmericanOptionTest::testFdAmericanGreeks() {
568	BOOST_TEST_MESSAGE("Testing finite-differences American option greeks...");
569	testFdGreeks<FdBlackScholesVanillaEngine>();
570	}
571
572	void AmericanOptionTest::testFdShoutGreeks() {
573	BOOST_TEST_MESSAGE("Testing finite-differences shout option greeks...");
574	testFdGreeks<FdBlackScholesShoutEngine>();
575	}
576
577	void AmericanOptionTest::testFDShoutNPV() {
578	BOOST_TEST_MESSAGE("Testing finite-differences shout option pricing...");
579
580	const auto dc = Actual365Fixed();
581	const auto today = Date(4, February, 2021);
582	Settings::instance().evaluationDate() = today;
583
584	const auto spot = Handle<Quote>(ext::make_shared<SimpleQuote>(args: 100.0));
585	const auto q = Handle<YieldTermStructure>(flatRate(forward: 0.03, dc));
586	const auto r = Handle<YieldTermStructure>(flatRate(forward: 0.06, dc));
587
588	const auto volTS = Handle<BlackVolTermStructure>(flatVol(volatility: 0.25, dc));
589	const auto process = ext::make_shared<BlackScholesMertonProcess>(
590	args: spot, args: q, args: r, args: volTS);
591
592	const auto maturityDate = today + Period(5, Years);
593
594	struct TestDescription { Real strike; Option::Type type; Real expected; };
595
596	const TestDescription testDescriptions[] = {
597	{.strike: 105, .type: Option::Put, .expected: 19.136},
598	{.strike: 105, .type: Option::Call, .expected: 28.211},
599	{.strike: 120, .type: Option::Put, .expected: 28.02},
600	{.strike: 80, .type: Option::Call, .expected: 40.785}
601	};
602
603	const auto engine = ext::make_shared<FdBlackScholesShoutEngine>(
604	args: process, args: 400, args: 200);
605
606	for (const TestDescription& desc: testDescriptions) {
607	const Real strike = desc.strike;
608	const Option::Type type = desc.type;
609
610	auto option = VanillaOption(
611	ext::make_shared<PlainVanillaPayoff>(args: type, args: strike),
612	ext::make_shared<AmericanExercise>(args: maturityDate));
613
614	option.setPricingEngine(engine);
615
616	const Real expected = desc.expected;
617	const Real tol = 2e-2;
618	const Real calculated = option.NPV();
619	const Real diff = std::fabs(x: calculated-expected);
620
621	if (diff > tol) {
622	BOOST_FAIL("failed to reproduce known shout option price for "
623	<< "\n strike: " << strike
624	<< "\n option type:" <<
625	((type == Option::Call)?"Call" : "Put")
626	<< "\n calculated: " << calculated
627	<< "\n expected: " << expected
628	<< "\n difference: " << diff
629	<< "\n tolerance: " << tol);
630	}
631	}
632	}
633
634	void AmericanOptionTest::testZeroVolFDShoutNPV() {
635	BOOST_TEST_MESSAGE("Testing zero volatility shout option pricing with discrete dividends...");
636
637	const auto dc = Actual365Fixed();
638	const auto today = Date(14, February, 2021);
639	Settings::instance().evaluationDate() = today;
640
641	const auto spot = Handle<Quote>(ext::make_shared<SimpleQuote>(args: 100.0));
642	const auto q = Handle<YieldTermStructure>(flatRate(forward: 0.03, dc));
643	const auto r = Handle<YieldTermStructure>(flatRate(forward: 0.07, dc));
644
645	const auto volTS = Handle<BlackVolTermStructure>(flatVol(volatility: 1e-6, dc));
646	const auto process = ext::make_shared<BlackScholesMertonProcess>(
647	args: spot, args: q, args: r, args: volTS);
648
649	const auto maturityDate = today + Period(1, Years);
650	const Date dividendDate = today + Period(3, Months);
651	const Real dividendAmount = 10.0;
652	auto dividends = DividendVector(dividendDates: { dividendDate },dividends: { dividendAmount });
653
654	VanillaOption option(
655	ext::make_shared<PlainVanillaPayoff>(args: Option::Put, args: 100.0),
656	ext::make_shared<AmericanExercise>(args: today, args: maturityDate)
657	);
658
659	option.setPricingEngine(
660	ext::make_shared<FdBlackScholesVanillaEngine>(
661	args: process, args&: dividends, args: 50, args: 50));
662
663	const Real americanNPV = option.NPV();
664
665	QL_DEPRECATED_DISABLE_WARNING
666	DividendVanillaOption divOption(
667	ext::make_shared<PlainVanillaPayoff>(args: Option::Put, args: 100.0),
668	ext::make_shared<AmericanExercise>(args: today, args: maturityDate),
669	std::vector<Date>{dividendDate},
670	std::vector<Real>{dividendAmount}
671	);
672	QL_DEPRECATED_ENABLE_WARNING
673
674	divOption.setPricingEngine(
675	ext::make_shared<FdBlackScholesShoutEngine>(args: process, args: 50, args: 50));
676
677	Real shoutNPV = divOption.NPV();
678	const DiscountFactor df = r->discount(d: maturityDate)/r->discount(d: dividendDate);
679
680	const Real tol = 1e-3;
681	Real diff = std::fabs(x: americanNPV - shoutNPV/df);
682
683	if (diff > tol) {
684	BOOST_FAIL("failed to reproduce American option NPV with "
685	"shout option pricing engine for "
686	<< "\n calculated: " << shoutNPV/df
687	<< "\n expected : " << americanNPV
688	<< "\n difference: " << diff
689	<< "\n tolerance: " << tol);
690	}
691
692	VanillaOption option2(
693	ext::make_shared<PlainVanillaPayoff>(args: Option::Put, args: 100.0),
694	ext::make_shared<AmericanExercise>(args: today, args: maturityDate)
695	);
696
697	option2.setPricingEngine(
698	ext::make_shared<FdBlackScholesShoutEngine>(args: process, args&: dividends, args: 50, args: 50));
699
700	shoutNPV = option2.NPV();
701	diff = std::fabs(x: americanNPV - shoutNPV/df);
702
703	if (diff > tol) {
704	BOOST_FAIL("failed to reproduce American option NPV with "
705	"shout option pricing engine for "
706	<< "\n calculated: " << shoutNPV/df
707	<< "\n expected : " << americanNPV
708	<< "\n difference: " << diff
709	<< "\n tolerance: " << tol);
710	}
711	}
712
713	void AmericanOptionTest::testLargeDividendShoutNPV() {
714	BOOST_TEST_MESSAGE("Testing zero strike shout option pricing with discrete dividends...");
715
716	const auto dc = Actual365Fixed();
717	const auto today = Date(21, February, 2021);
718	Settings::instance().evaluationDate() = today;
719
720	const Real s0 = 100.0;
721	const Volatility vol = 0.25;
722
723	const auto q = Handle<YieldTermStructure>(flatRate(forward: 0.00, dc));
724	const auto r = Handle<YieldTermStructure>(flatRate(forward: 0.00, dc));
725	const auto vTS = Handle<BlackVolTermStructure>(flatVol(volatility: vol, dc));
726
727	const auto process = ext::make_shared<BlackScholesMertonProcess>(
728	args: Handle<Quote>(ext::make_shared<SimpleQuote>(args: s0)), args: q, args: r, args: vTS);
729
730	const auto maturityDate = today + Period(6, Months);
731	const Date dividendDate = today + Period(3, Months);
732	const Real divAmount = 30.0;
733	auto dividends = DividendVector(dividendDates: { dividendDate }, dividends: { divAmount });
734
735	const Real strike = 80.0;
736	VanillaOption option(
737	ext::make_shared<PlainVanillaPayoff>(args: Option::Call, args: strike),
738	ext::make_shared<AmericanExercise>(args: today, args: maturityDate)
739	);
740
741	option.setPricingEngine(
742	ext::make_shared<FdBlackScholesShoutEngine>(args: process, args&: dividends, args: 100, args: 400));
743
744	Real calculated = option.NPV();
745
746	VanillaOption ref_option(
747	ext::make_shared<PlainVanillaPayoff>(args: Option::Call, args: strike),
748	ext::make_shared<AmericanExercise>(args: today, args: dividendDate)
749	);
750
751	ref_option.setPricingEngine(
752	ext::make_shared<FdBlackScholesShoutEngine>(args: process, args: 100, args: 400));
753
754	const Real expected = ref_option.NPV()
755	* r->discount(d: maturityDate) / r->discount(d: dividendDate);
756
757	const Real tol = 5e-2;
758	Real diff = std::fabs(x: expected - calculated);
759
760	if (diff > tol) {
761	BOOST_FAIL("failed to reproduce American option NPV with "
762	"shout option pricing engine for "
763	<< "\n calculated: " << calculated
764	<< "\n expected : " << expected
765	<< "\n difference: " << diff
766	<< "\n tolerance: " << tol);
767	}
768
769	QL_DEPRECATED_DISABLE_WARNING
770	DividendVanillaOption divOption(
771	ext::make_shared<PlainVanillaPayoff>(args: Option::Call, args: strike),
772	ext::make_shared<AmericanExercise>(args: today, args: maturityDate),
773	std::vector<Date>{dividendDate},
774	std::vector<Real>{divAmount}
775	);
776	QL_DEPRECATED_ENABLE_WARNING
777
778	divOption.setPricingEngine(
779	ext::make_shared<FdBlackScholesShoutEngine>(args: process, args: 100, args: 400));
780
781	calculated = divOption.NPV();
782
783	diff = std::fabs(x: expected - calculated);
784
785	if (diff > tol) {
786	BOOST_FAIL("failed to reproduce American option NPV with "
787	"shout option pricing engine for "
788	<< "\n calculated: " << calculated
789	<< "\n expected : " << expected
790	<< "\n difference: " << diff
791	<< "\n tolerance: " << tol);
792	}
793	}
794
795	void AmericanOptionTest::testEscrowedVsSpotAmericanOption() {
796	BOOST_TEST_MESSAGE("Testing escrowed vs spot dividend model for American options...");
797
798	const auto dc = Actual360();
799	const auto today = Date(27, February, 2021);
800	Settings::instance().evaluationDate() = today;
801
802	ext::shared_ptr<SimpleQuote> vol(new SimpleQuote(0.3));
803
804	const auto process = ext::make_shared<BlackScholesMertonProcess>(
805	args: Handle<Quote>(ext::make_shared<SimpleQuote>(args: 100)),
806	args: Handle<YieldTermStructure>(flatRate(forward: 0.08, dc)),
807	args: Handle<YieldTermStructure>(flatRate(forward: 0.04, dc)),
808	args: Handle<BlackVolTermStructure>(flatVol(volatility: vol, dc))
809	);
810
811	const auto maturityDate = today + Period(12, Months);
812	std::vector<Date> dividendDates = { today + Period(10, Months) };
813	std::vector<Real> dividendAmounts = { 10.0 };
814	auto dividends = DividendVector(dividendDates, dividends: dividendAmounts);
815
816	const Real strike = 100.0;
817	VanillaOption option(
818	ext::make_shared<PlainVanillaPayoff>(args: Option::Call, args: strike),
819	ext::make_shared<AmericanExercise>(args: today, args: maturityDate)
820	);
821
822	option.setPricingEngine(
823	ext::make_shared<FdBlackScholesVanillaEngine>(args: process, args&: dividends, args: 100, args: 400));
824
825	const Real spotNpv = option.NPV();
826	const Real spotDelta = option.delta();
827
828	vol->setValue(100/90.*0.3);
829
830	option.setPricingEngine(
831	MakeFdBlackScholesVanillaEngine(process)
832	.withTGrid(tGrid: 100)
833	.withXGrid(xGrid: 400)
834	.withCashDividends(dividendDates, dividendAmounts)
835	.withCashDividendModel(cashDividendModel: FdBlackScholesVanillaEngine::Escrowed)
836	);
837
838	const Real escrowedNpv = option.NPV();
839	const Real escrowedDelta = option.delta();
840
841	const Real diffNpv = std::abs(x: escrowedNpv - spotNpv);
842	const Real tol = 1e-2;
843
844	if (diffNpv > tol) {
845	BOOST_FAIL("failed to compare American option NPV with "
846	"escrowed and spot dividend model "
847	<< "\n escrowed div: " << escrowedNpv
848	<< "\n spot div : " << spotNpv
849	<< "\n difference: " << diffNpv
850	<< "\n tolerance: " << tol);
851	}
852
853
854	const Real diffDelta = std::abs(x: escrowedDelta - spotDelta);
855
856	if (diffDelta > tol) {
857	BOOST_FAIL("failed to compare American option Delta with "
858	"escrowed and spot dividend model "
859	<< "\n escrowed div: " << escrowedDelta
860	<< "\n spot div : " << spotDelta
861	<< "\n difference: " << diffDelta
862	<< "\n tolerance: " << tol);
863	}
864	}
865
866
867	void AmericanOptionTest::testTodayIsDividendDate() {
868	BOOST_TEST_MESSAGE("Testing escrowed vs spot dividend model on dividend dates for American options...");
869
870	const auto dc = Actual360();
871	const auto today = Date(27, February, 2021);
872	Settings::instance().evaluationDate() = today;
873
874	ext::shared_ptr<SimpleQuote> vol(new SimpleQuote(0.3));
875
876	const auto process = ext::make_shared<BlackScholesMertonProcess>(
877	args: Handle<Quote>(ext::make_shared<SimpleQuote>(args: 100)),
878	args: Handle<YieldTermStructure>(flatRate(forward: 0.05, dc)),
879	args: Handle<YieldTermStructure>(flatRate(forward: 0.07, dc)),
880	args: Handle<BlackVolTermStructure>(flatVol(volatility: vol, dc))
881	);
882
883	const auto maturityDate = today + Period(12, Months);
884	std::vector<Date> dividendDates = { today, today + Period(11, Months) };
885	std::vector<Real> dividendAmounts = { 5.0, 5.0 };
886
887	ext::shared_ptr<PricingEngine> spotEngine =
888	MakeFdBlackScholesVanillaEngine(process)
889	.withTGrid(tGrid: 100)
890	.withXGrid(xGrid: 400)
891	.withCashDividends(dividendDates, dividendAmounts)
892	.withCashDividendModel(cashDividendModel: FdBlackScholesVanillaEngine::Spot);
893
894	ext::shared_ptr<PricingEngine> escrowedEngine =
895	MakeFdBlackScholesVanillaEngine(process)
896	.withTGrid(tGrid: 100)
897	.withXGrid(xGrid: 400)
898	.withCashDividends(dividendDates, dividendAmounts)
899	.withCashDividendModel(cashDividendModel: FdBlackScholesVanillaEngine::Escrowed);
900
901	const Real strike = 90.0;
902	VanillaOption option(
903	ext::make_shared<PlainVanillaPayoff>(args: Option::Put, args: strike),
904	ext::make_shared<AmericanExercise>(args: today, args: maturityDate)
905	);
906
907	option.setPricingEngine(spotEngine);
908
909	Real spotNpv = option.NPV();
910	const Real spotDelta = option.delta();
911	BOOST_CHECK_THROW(option.theta(), QuantLib::Error);
912
913	vol->setValue(100/95.*0.3);
914
915	option.setPricingEngine(escrowedEngine);
916
917	Real escrowedNpv = option.NPV();
918	const Real escrowedDelta = option.delta();
919	BOOST_CHECK_THROW(option.theta(), QuantLib::Error);
920
921	Real diffNpv = std::abs(x: escrowedNpv - spotNpv);
922	Real tol = 5e-2;
923
924	if (diffNpv > tol) {
925	BOOST_FAIL("failed to compare American option NPV with "
926	"escrowed and spot dividend model "
927	<< "\n escrowed div: " << escrowedNpv
928	<< "\n spot div : " << spotNpv
929	<< "\n difference: " << diffNpv
930	<< "\n tolerance: " << tol);
931	}
932
933	const Real diffDelta = std::abs(x: escrowedDelta - spotDelta);
934
935	tol = 1e-3;
936	if (diffDelta > tol) {
937	BOOST_FAIL("failed to compare American option Delta with "
938	"escrowed and spot dividend model "
939	<< "\n escrowed div: " << escrowedDelta
940	<< "\n spot div : " << spotDelta
941	<< "\n difference: " << diffDelta
942	<< "\n tolerance: " << tol);
943	}
944
945	dividendDates[0] = today + 1;
946
947	spotEngine =
948	MakeFdBlackScholesVanillaEngine(process)
949	.withTGrid(tGrid: 100)
950	.withXGrid(xGrid: 400)
951	.withCashDividends(dividendDates, dividendAmounts)
952	.withCashDividendModel(cashDividendModel: FdBlackScholesVanillaEngine::Spot);
953
954	escrowedEngine =
955	MakeFdBlackScholesVanillaEngine(process)
956	.withTGrid(tGrid: 100)
957	.withXGrid(xGrid: 400)
958	.withCashDividends(dividendDates, dividendAmounts)
959	.withCashDividendModel(cashDividendModel: FdBlackScholesVanillaEngine::Escrowed);
960
961	vol->setValue(0.3);
962
963	option.setPricingEngine(spotEngine);
964	spotNpv = option.NPV();
965
966	vol->setValue(100/95.0*0.3);
967	option.setPricingEngine(escrowedEngine);
968
969	escrowedNpv = option.NPV();
970	BOOST_CHECK_NO_THROW(option.theta());
971
972	diffNpv = std::abs(x: escrowedNpv - spotNpv);
973	tol = 5e-2;
974
975	if (diffNpv > tol) {
976	BOOST_FAIL("failed to compare American option NPV with "
977	"escrowed and spot dividend model "
978	<< "\n escrowed div: " << escrowedNpv
979	<< "\n spot div : " << spotNpv
980	<< "\n difference: " << diffNpv
981	<< "\n tolerance: " << tol);
982	}
983	}
984
985
986	void AmericanOptionTest::testCallPutParity() {
987	BOOST_TEST_MESSAGE("Testing call/put parity for American options...");
988
989	// R.L. McDonald, M.D. Schroder: A parity result for American option
990
991	const DayCounter dc = Actual365Fixed();
992	const Date today = Date(8, April, 2022);
993	Settings::instance().evaluationDate() = today;
994
995	struct OptionSpec {
996	Real spot;
997	Real strike;
998	Size maturityInDays;
999	Real volatility;
1000	Real r;
1001	Real q;
1002	};
1003
1004	auto buildStochProcess = [&dc](const OptionSpec& testCase) {
1005	return ext::make_shared<BlackScholesMertonProcess>(
1006	args: Handle<Quote>(ext::make_shared<SimpleQuote>(args: testCase.spot)),
1007	args: Handle<YieldTermStructure>(flatRate(forward: testCase.q, dc)),
1008	args: Handle<YieldTermStructure>(flatRate(forward: testCase.r, dc)),
1009	args: Handle<BlackVolTermStructure>(flatVol(volatility: testCase.volatility, dc))
1010	);
1011	};
1012	const OptionSpec testCaseSpecs[] = {
1013	{.spot: 100.0, .strike: 100.0, .maturityInDays: 365, .volatility: 0.5, .r: 0.15, .q: 0.02},
1014	{.spot: 100.0, .strike: 90.0, .maturityInDays: 365, .volatility: 0.5, .r: 0.15, .q: 0.02},
1015	{.spot: 100.0, .strike: 125.0, .maturityInDays: 730, .volatility: 0.4, .r: 0.15, .q: 0.05},
1016	{.spot: 100.0, .strike: 125.0, .maturityInDays: 730, .volatility: 0.4, .r: 0.06, .q: 0.05}
1017	};
1018
1019	const Size xGrid = 400;
1020	const Size timeStepsPerYear=50;
1021
1022	for (const auto& testCaseSpec: testCaseSpecs) {
1023	const auto maturityDate =
1024	today + Period(testCaseSpec.maturityInDays, Days);
1025	const Time maturityTime = dc.yearFraction(d1: today, d2: maturityDate);
1026	const Size tGrid = Size(maturityTime * timeStepsPerYear);
1027
1028	const auto exercise =
1029	ext::make_shared<AmericanExercise>(args: today, args: maturityDate);
1030
1031	VanillaOption putOption(
1032	ext::make_shared<PlainVanillaPayoff>(
1033	args: Option::Put, args: testCaseSpec.strike),
1034	exercise
1035	);
1036	putOption.setPricingEngine(
1037	ext::make_shared<FdBlackScholesVanillaEngine>(
1038	args: buildStochProcess(testCaseSpec), args: tGrid, args: xGrid)
1039	);
1040	const Real putNpv = putOption.NPV();
1041
1042	OptionSpec callOptionSpec = {
1043	.spot: testCaseSpec.strike,
1044	.strike: testCaseSpec.spot,
1045	.maturityInDays: testCaseSpec.maturityInDays,
1046	.volatility: testCaseSpec.volatility,
1047	.r: testCaseSpec.q,
1048	.q: testCaseSpec.r
1049	};
1050	VanillaOption callOption(
1051	ext::make_shared<PlainVanillaPayoff>(
1052	args: Option::Call, args&: callOptionSpec.strike),
1053	exercise
1054	);
1055	callOption.setPricingEngine(
1056	ext::make_shared<FdBlackScholesVanillaEngine>(
1057	args: buildStochProcess(callOptionSpec), args: tGrid, args: xGrid)
1058	);
1059	const Real callNpv = callOption.NPV();
1060
1061	const Real diff = std::fabs(x: putNpv -callNpv);
1062	const Real tol = 0.001;
1063
1064	if (diff > tol) {
1065	BOOST_FAIL("failed to reproduce American call/put parity"
1066	<< "\n Put NPV : " << putNpv
1067	<< "\n Call NPV : " << callNpv
1068	<< "\n difference: " << diff
1069	<< "\n tolerance : " << tol);
1070	}
1071	}
1072	}
1073
1074	void AmericanOptionTest::testQdPlusBoundaryValues() {
1075	BOOST_TEST_MESSAGE("Testing QD+ boundary approximation...");
1076
1077	const DayCounter dc = Actual365Fixed();
1078	const Real S = 100;
1079	const Real K = 120;
1080	const Rate r = 0.1;
1081	const Rate q = 0.03;
1082	const Volatility sigma = 0.25;
1083	const Time maturity = 5.0;
1084
1085	const QdPlusAmericanEngine qrPlusEngine(
1086	ext::shared_ptr<GeneralizedBlackScholesProcess>(), 10);
1087
1088	std::vector<std::pair<Real, Real> > testCaseSpecs = {
1089	{4.9, 87.76960949965387},
1090	{4.0, 88.39053003614612},
1091	{2.5, 90.14327315762256},
1092	{1.0, 94.49793803095984},
1093	{0.1, 106.2588964442338}
1094	};
1095
1096	for (const auto& testCaseSpec: testCaseSpecs) {
1097	const auto calculated
1098	= qrPlusEngine.putExerciseBoundaryAtTau(
1099	S, K, r, q, vol: sigma, T: maturity, tau: testCaseSpec.first);
1100
1101	const Real boundary = calculated.second;
1102	const Size nrEvaluations = calculated.first;
1103
1104	const Real expected = testCaseSpec.second;
1105
1106	const Real diff = std::fabs(x: boundary - expected);
1107	const Real tol = 1e-12;
1108
1109	if (diff > tol) {
1110	BOOST_FAIL("failed to reproduce QR+ boundary approximation"
1111	<< "\n calculated: " << boundary
1112	<< "\n expected: " << expected
1113	<< "\n difference: " << diff
1114	<< "\n tolerance : " << tol);
1115	}
1116
1117	if (nrEvaluations > 10) {
1118	BOOST_FAIL("failed to reproduce rate of convergence"
1119	<< "\n evaluations: " << nrEvaluations
1120	<< "\n max eval : " << 10);
1121	}
1122	}
1123	}
1124
1125	void AmericanOptionTest::testQdPlusBoundaryConvergence() {
1126	BOOST_TEST_MESSAGE("Testing QD+ boundary convergence...");
1127
1128	const DayCounter dc = Actual365Fixed();
1129	const Real S = 100;
1130	const Volatility sigma = 0.25;
1131	const Time maturity = 10.0;
1132
1133	struct TestCaseSpec {
1134	Real r, q, strike;
1135	Size maxEvaluations;
1136	};
1137
1138	TestCaseSpec testCases[] = {
1139	{ .r: 0.10, .q: 0.03, .strike: 120, .maxEvaluations: 2000 },
1140	{ .r: 0.0001, .q: 0.03, .strike: 120, .maxEvaluations: 2000 },
1141	{ .r: 0.0001, .q: 0.000002, .strike: 120, .maxEvaluations: 2000 },
1142	{ .r: 0.01, .q: 0.75, .strike: 120, .maxEvaluations: 2000 },
1143	{ .r: 0.03, .q: 0.0, .strike: 30, .maxEvaluations: 2000 },
1144	{ .r: 0.03, .q: 0.0, .strike: 1e7, .maxEvaluations: 2500 },
1145	{ .r: 0.075, .q: 0.0, .strike: 1e-8, .maxEvaluations: 2000 }
1146	};
1147
1148	const std::vector<std::pair<QdPlusAmericanEngine::SolverType, std::string> >
1149	solverTypes{
1150	{ QdPlusAmericanEngine::Brent, "Brent" },
1151	{ QdPlusAmericanEngine::Newton, "Newton" },
1152	{ QdPlusAmericanEngine::Ridder, "Ridder" },
1153	{ QdPlusAmericanEngine::Halley, "Halley" },
1154	{ QdPlusAmericanEngine::SuperHalley, "SuperHalley" }
1155	};
1156
1157	for (const auto& testCase: testCases) {
1158	for (const auto& solverType : solverTypes) {
1159	const QdPlusAmericanEngine qrPlusEngine(
1160	ext::shared_ptr<GeneralizedBlackScholesProcess>(),
1161	Null<Size>(), solverType.first, 1e-8);
1162
1163	Size nrEvaluations = 0;
1164
1165	for (Real t=0.0; t < maturity; t+=0.1) {
1166	const auto calculated = qrPlusEngine.putExerciseBoundaryAtTau(
1167	S, K: testCase.strike, r: testCase.r,
1168	q: testCase.q, vol: sigma, T: maturity, tau: t);
1169	nrEvaluations += calculated.first;
1170	}
1171
1172	const Size maxEvaluations =
1173	( solverType.first == QdPlusAmericanEngine::Halley
1174	|| solverType.first == QdPlusAmericanEngine::SuperHalley)
1175	? 750 : testCase.maxEvaluations;
1176
1177	if (nrEvaluations > maxEvaluations) {
1178	BOOST_FAIL("QR+ boundary approximation failed to converge "
1179	<< "\n evaluations: " << nrEvaluations
1180	<< "\n max eval: " << maxEvaluations
1181	<< "\n Solver: " << solverType.second
1182	<< "\n r : " << testCase.r
1183	<< "\n q : " << testCase.q
1184	<< "\n K : " << testCase.strike);
1185	}
1186	}
1187	}
1188	}
1189
1190	void AmericanOptionTest::testQdAmericanEngines() {
1191	BOOST_TEST_MESSAGE("Testing QD+ American option pricing...");
1192
1193	const DayCounter dc = Actual365Fixed();
1194	const Date today = Date(1, June, 2022);
1195	Settings::instance().evaluationDate() = today;
1196
1197	struct OptionSpec {
1198	Option::Type optionType;
1199	Real spot;
1200	Real strike;
1201	Size maturityInDays;
1202	Real volatility;
1203	Real r;
1204	Real q;
1205	Real expectedValue;
1206	Real precision;
1207	};
1208
1209	// high precision edge cases
1210	const OptionSpec edgeTestCases[] = {
1211
1212	// standard put option
1213	{.optionType: Option::Put, .spot: 100.0, .strike: 120.0, .maturityInDays: 3650, .volatility: 0.25, .r: 0.10, .q: 0.03, .expectedValue: 22.97383256003585, .precision: 1e-8},
1214	// call-put parity on standard option
1215	{.optionType: Option::Call, .spot: 120.0, .strike: 100.0, .maturityInDays: 3650, .volatility: 0.25, .r: 0.03, .q: 0.10, .expectedValue: 22.97383256003585, .precision: 1e-8},
1216
1217	// zero strike put
1218	{.optionType: Option::Put, .spot: 100.0, .strike: 0.0, .maturityInDays: 365, .volatility: 0.25, .r: 0.02, .q: 0.02, .expectedValue: 0.0, .precision: 1e-14},
1219	{.optionType: Option::Put, .spot: 100.0, .strike: 1e-8, .maturityInDays: 365, .volatility: 0.25, .r: 0.02, .q: 0.02, .expectedValue: 0.0, .precision: 1e-14},
1220
1221	// zero strike call
1222	{.optionType: Option::Call, .spot: 100.0, .strike: 0.0, .maturityInDays: 365, .volatility: 0.25, .r: 0.05, .q: 0.01, .expectedValue: 100.0, .precision: 1e-11},
1223	{.optionType: Option::Call, .spot: 100.0, .strike: 1e-7, .maturityInDays: 365, .volatility: 0.25, .r: 0.05, .q: 0.01, .expectedValue: 100.0-1e-7, .precision: 1e-9},
1224
1225	// zero vol call
1226	{.optionType: Option::Call, .spot: 100.0, .strike: 50.0, .maturityInDays: 365, .volatility: 0.0, .r: 0.05, .q: 0.01, .expectedValue: 51.4435121498811085, .precision: 1e-10},
1227	{.optionType: Option::Call, .spot: 100.0, .strike: 50.0, .maturityInDays: 365, .volatility: 1e-8, .r: 0.05, .q: 0.01, .expectedValue: 51.4435121498811156, .precision: 1e-8},
1228
1229	// zero vol put 1
1230	{.optionType: Option::Put, .spot: 100.0, .strike: 120.0, .maturityInDays: 4*3650, .volatility: 1e-6, .r: 0.01, .q: 0.50, .expectedValue: 108.980920365700442, .precision: 1e-4},
1231	{.optionType: Option::Put, .spot: 100.0, .strike: 120.0, .maturityInDays: 4*3650, .volatility: 0.0, .r: 0.01, .q: 0.50, .expectedValue: 108.980904561184602, .precision: 1e-10},
1232
1233	// zero vol put 2
1234	{.optionType: Option::Put, .spot: 100.0, .strike: 120.0, .maturityInDays: 365, .volatility: 1e-7, .r: 0.05, .q: 0.01, .expectedValue: 20.0, .precision: 1e-9},
1235	{.optionType: Option::Put, .spot: 100.0, .strike: 120.0, .maturityInDays: 365, .volatility: 0.0, .r: 0.05, .q: 0.01, .expectedValue: 20.0, .precision: 1e-12},
1236
1237	// zero vol put 3
1238	{.optionType: Option::Put, .spot: 100.0, .strike: 120.0, .maturityInDays: 365, .volatility: 1e-7, .r: 0.00, .q: 0.05, .expectedValue: 24.8770575499286082, .precision: 1e-8},
1239	{.optionType: Option::Put, .spot: 100.0, .strike: 120.0, .maturityInDays: 365, .volatility: 0.0, .r: 0.00, .q: 0.05, .expectedValue: 24.8770575499286082, .precision: 1e-10},
1240
1241	// zero spot put
1242	{.optionType: Option::Put, .spot: 1e-6, .strike: 120.0, .maturityInDays: 365, .volatility: 0.25, .r: -0.075, .q: 0.05, .expectedValue: 129.346097154926355, .precision: 1e-9},
1243	{.optionType: Option::Put, .spot: 0.0, .strike: 120.0, .maturityInDays: 365, .volatility: 0.25, .r: -0.075, .q: 0.05, .expectedValue: 129.346098106155779, .precision: 1e-10},
1244
1245	// zero spot call
1246	{.optionType: Option::Call, .spot: 1e-6, .strike: 120.0, .maturityInDays: 365, .volatility: 0.25, .r: 0.075, .q: 0.05, .expectedValue: 0.0, .precision: 1e-14},
1247	{.optionType: Option::Call, .spot: 0.0, .strike: 120.0, .maturityInDays: 365, .volatility: 0.25, .r: 0.075, .q: 0.05, .expectedValue: 0.0, .precision: 1e-14},
1248
1249	// put option with one day left
1250	{.optionType: Option::Put, .spot: 100.0, .strike: 120.0, .maturityInDays: 1, .volatility: 0.25, .r: 0.05, .q: 0.0, .expectedValue: 20.0, .precision: 1e-10},
1251
1252	// put option at maturity
1253	{.optionType: Option::Put, .spot: 100.0, .strike: 120.0, .maturityInDays: 0, .volatility: 0.25, .r: 0.05, .q: 0.0, .expectedValue: 0.0, .precision: 1e-14},
1254
1255	// zero everything
1256	{.optionType: Option::Put, .spot: 0.0, .strike: 0.0, .maturityInDays: 365, .volatility: 0.0, .r: 0.0, .q: 0.0, .expectedValue: 0.0, .precision: 1e-14},
1257
1258	// zero strike call with zero vol
1259	{.optionType: Option::Call, .spot: 100.0, .strike: 1e-7, .maturityInDays: 365, .volatility: 1e-8, .r: 0.05, .q: 0.025, .expectedValue: 100.0-1e-7, .precision: 1e-8},
1260	{.optionType: Option::Call, .spot: 100.0, .strike: 0.0, .maturityInDays: 365, .volatility: 1e-8, .r: 0.05, .q: 0.025, .expectedValue: 100.0, .precision: 1e-8},
1261	{.optionType: Option::Call, .spot: 100.0, .strike: 1e-7, .maturityInDays: 365, .volatility: 0.0, .r: 0.05, .q: 0.025, .expectedValue: 100.0-1e-7, .precision: 1e-8},
1262	{.optionType: Option::Call, .spot: 100.0, .strike: 0.0, .maturityInDays: 365, .volatility: 0.0, .r: 0.05, .q: 0.025, .expectedValue: 100.0, .precision: 1e-8},
1263
1264	// zero spot call with zero vol
1265	{.optionType: Option::Call, .spot: 1e-8, .strike: 100, .maturityInDays: 365, .volatility: 1e-8, .r: 0.05, .q: 0.025, .expectedValue: 0.0, .precision: 1e-10},
1266	{.optionType: Option::Call, .spot: 0.0, .strike: 100, .maturityInDays: 365, .volatility: 0.0, .r: 0.05, .q: 0.025, .expectedValue: 0.0, .precision: 1e-14},
1267
1268	// zero interest rate call
1269	{.optionType: Option::Call, .spot: 100, .strike: 100, .maturityInDays: 365, .volatility: 0.25, .r: 0.0, .q: 0.025, .expectedValue: 8.871505915120776, .precision: 1e-8},
1270
1271	// zero dividend rate call
1272	{.optionType: Option::Call, .spot: 100, .strike: 100, .maturityInDays: 365, .volatility: 0.25, .r: 0.05, .q: 0.0, .expectedValue: 12.3359989303687243, .precision: 1e-8},
1273
1274	// extreme spot call
1275	{.optionType: Option::Call, .spot: 1e10, .strike: 100, .maturityInDays: 365, .volatility: 0.25, .r: 0.01, .q: 0.05, .expectedValue: 1e10-100.0, .precision: -1},
1276
1277	// extreme strike call
1278	{.optionType: Option::Call, .spot: 100, .strike: 1e10, .maturityInDays: 365, .volatility: 0.25, .r: 0.01, .q: 0.05, .expectedValue: 0.0, .precision: 1e-14},
1279
1280	// extreme vol call
1281	{.optionType: Option::Call, .spot: 100, .strike: 100, .maturityInDays: 365, .volatility: 100.0, .r: 0.01, .q: 0.05, .expectedValue: 99.9874942266127, .precision: 1e-8},
1282
1283	// extreme dividend yield call
1284	{.optionType: Option::Call, .spot: 100, .strike: 100, .maturityInDays: 365, .volatility: 0.25, .r: 0.10, .q: 10.0, .expectedValue: 0.1159627202107989, .precision: 1e-8},
1285
1286	// extreme maturity call
1287	{.optionType: Option::Call, .spot: 100, .strike: 100, .maturityInDays: 170*365, .volatility: 0.25, .r: 0.01, .q: 0.002, .expectedValue: 80.37468392429741, .precision: 1e-8}
1288	};
1289
1290	// random test cases
1291	const double pde_values[] = {
1292	581.46895,113.78442,581.44547,1408.579,49.19448,1060.27367,
1293	834.83366,176.48305,120.38008,307.11264,602.7006,233.80171,
1294	204.74596,0.30987,0,0,5.36215,0.01711,0,84.51193,0.67131,
1295	0.06414,152.67188,54.75257,90.31861,168.50289,18.38926,0,
1296	282.4995,0,0.08428,12.30929,42.26359,139.87748,0.28724,0.00421,
1297	0,0.00206,0,658.60427,140.51139,23.17387,0.35612,0,909.14828,
1298	0,0.11549,5.46749,144.25428,2576.6754,562.16484,0,122.725,
1299	383.48463,278.7447,3.52566,82.34348,81.06139,0,10.42824,
1300	4.95917,25.28602,31.38869,3.53697,0,0.012,0,0.4263,162.16184,
1301	0.4618,97.714,283.03442,0.38176,70.25367,134.94142,2.19293,
1302	226.4746,76.74309,46.03123,15.76214,0.01666,1806.26208,0,
1303	103.93726,6.82956,337.81301,0.64236,677.63248,25.01763,
1304	443.79052,1793.78327,118.6293,185.79849,11.59313,679.01736,
1305	17.99005,403.57554,1.67418,0,0.03795,3326.09089,71.1996,
1306	0,485.10353,0,1681.25166,0,43.15432,0.75825,0.05895,34.71493,
1307	0.00015,5.58671,115.98793,37.7713,399.24494,0.00766,445.42207,
1308	152.65397,0,47.05874,0.96921,14.21875,257.84754,109.62533,
1309	2553.99295,138.46663,192.33614,81.41877,18.21403,113.926,
1310	27.28409,174.77093,42.70527,0.90326,0,967.9901,616.0143,
1311	253.56442,0.00397,2493.82098,9.29406,11.00023,0,0,234.12481,
1312	0,72.46356,0,9.00932,48.67934,29.42756,13.4271,0,0,0,0,20.71417,
1313	48.57474,2.26452,0,109.0243,0,21.26801,1.21164,0,86.25232,
1314	36.00437,4.53844,7.40503,313.53602,379.76105,165.84347,77.19665,
1315	9.02466,0.10634,214.84982,6.13387,133.44645,303.25953,0,
1316	134.26724,246.89804,0,123.32975,32.83429,9.56819,7.42582,0,
1317	73.82832,196.84831,0.00001,72.70391,2173.8649,123.00513,
1318	153.83539,21.63003,209.84752,30.12425,0,197.6502,0,164.02863,
1319	7.65143,56.57631,2392.70018,0,0,34.23457,171.08459,0.49387,
1320	31.13395,237.68801,0.01262,0,0,0,0,41.56635,0,8.41535,55.01775,
1321	310.50094,0,14.85456,174.34018,7.19772,0.00001,0,91.70874,
1322	0.00001,17.51724,0.00587,0,532.24902,2.05553,36.80843,0,
1323	33.39288,0.00006,0.04439,1.3434,0,0.41816,926.37642,0,247.61559,
1324	151.98965,0.35243,4.33198,23294.47744,0.00791,12.51996,53.47727,
1325	167.95572,0.0062,6.8482,0,347.83408,852.85742,558.21422,0,
1326	53.89293,78.61011,187.3978,9.18927,0.00553,113.48101,1467.30556,
1327	74.82251,94.84476,0,101.3649,59.27007,0,773.81251,0,542.7889,0,
1328	68.96209,96.0435,0.00004,0.10738,0.00187,324.97758,245.68455,
1329	30.52818,129.84472,0,46.86288,368.41675,139.29763,4.4393,16.29594,
1330	25.7554,64.02621,89.41363,0.62751,219.65237,0.26039,0,12.02172,
1331	101.97733,69.37456,45.81122,1263.33603,164.31607,15.88788,0,48.77797,
1332	0.13233,147.16808,10.31217,7.50634,7.48611,177.95409,225.77562,3.56947,
1333	0.02531,4.88869,8.76632,0,0,0.02214,305.08468,44.52185,182.17332,
1334	538.31458,0,46.97229,0,31.94202,410.43038,0,70.35432,15.58346,74.14177,
1335	953.67663,11.79128,59.83061,0,37.86557,1184.22731,2411.37823,0,0,0,0,
1336	49.3179,236.38654,21.36225,218.048,517.57006,0,0,12.52933,256.71967,
1337	0.00025,1.47981,158.19166,0,1923.70709,4.94441,1199.81196,45.92353,
1338	85.73255,14.91338,88.81459,21.42459,3456.9466,31.97838,233.26863,
1339	49.34801,2684.07758,0,0,32.24149,0,111.79552,0.00506,8.77602,0,
1340	406.54213,0.32974,365.53998,1.49714,19.65603,37.33877,205.06928,
1341	0.01805,589.23478,9.58273,0.02946,286.48706,463.34512,528.21392,0,
1342	47.71294,21.0864,114.81771,80.489,21.30905,41.95873,19.03598,156.09295,
1343	0,73.6509,0,0,168.17576,0,32.71243,36.75044,177.64583,0.05618,
1344	156.38616,1370.4754,24.5976,59.83173,0,354.93074,34.96889,0.00532,
1345	16.95287,1259.72993,241.05777,18.9778,0.57635,43.98093,25.2678,
1346	369.39896,0.31549,0,31.95512,101.60559,11.22079,970.16273,0,0,
1347	1.55445,0,18.6067,0,1124.20117,52.67762,10.38273,0,10.22588,251.27813,
1348	0,431.82244,0,1.31252,0,84.72154,100.98411,160.95557,129.51372,
1349	0.00026,103.81663,421.64767,0.00031,0,104.48529,162.59225,0,
1350	1504.0869,88.11253,4.14052,0.07195,203.78754,0.00002,42.5395,0,
1351	17.05087,26.89157,64.64923,0,390.87453,124.55406,0.01018,94.23963};
1352
1353	std::vector<OptionSpec> testCaseSpecs;
1354	testCaseSpecs.reserve(LENGTH(pde_values) + LENGTH(edgeTestCases));
1355
1356	PseudoRandom::rng_type rng(PseudoRandom::urng_type(12345UL));
1357
1358	for (double pde_value : pde_values) {
1359	const Option::Type optionType
1360	= (rng.next().value > 0)? Option::Call : Option::Put;
1361	const Real spot = 100*std::exp(x: 1.5*rng.next().value);
1362	const Real strike = 100*std::exp(x: 1.5*rng.next().value);
1363	const Size maturityInDays = Size(1 + 365*std::exp(x: 2*rng.next().value));
1364	const Volatility vol = 0.5*std::exp(x: rng.next().value);
1365	const Rate r = 0.10*std::exp(x: rng.next().value);
1366	const Rate q = 0.10*std::exp(x: rng.next().value);
1367
1368	const OptionSpec spec = {.optionType: optionType, .spot: spot, .strike: strike, .maturityInDays: maturityInDays, .volatility: vol, .r: r,
1369	.q: q, .expectedValue: pde_value, .precision: -1};
1370
1371	testCaseSpecs.push_back(x: spec);
1372	}
1373
1374	testCaseSpecs.insert(
1375	position: testCaseSpecs.end(),first: std::begin(arr: edgeTestCases), last: std::end(arr: edgeTestCases));
1376
1377	const auto spot = ext::make_shared<SimpleQuote>(args: 1.0);
1378	const auto rRate = ext::make_shared<SimpleQuote>(args: 0.0);
1379	const auto qRate = ext::make_shared<SimpleQuote>(args: 0.0);
1380	const auto vol = ext::make_shared<SimpleQuote>(args: 0.0);
1381
1382	const auto bsProcess = ext::make_shared<BlackScholesMertonProcess>(
1383	args: Handle<Quote>(spot),
1384	args: Handle<YieldTermStructure>(flatRate(today, forward: qRate, dc)),
1385	args: Handle<YieldTermStructure>(flatRate(today, forward: rRate, dc)),
1386	args: Handle<BlackVolTermStructure>(flatVol(today, volatility: vol, dc))
1387	);
1388
1389	const auto qrPlusAmericanEngine =
1390	ext::make_shared<QdPlusAmericanEngine>(
1391	args: bsProcess, args: 8, args: QdPlusAmericanEngine::Halley, args: 1e-10
1392	);
1393
1394	for (const auto& testCaseSpec: testCaseSpecs) {
1395	const Date maturityDate =
1396	today + Period(testCaseSpec.maturityInDays, Days);
1397
1398	spot->setValue(testCaseSpec.spot);
1399	rRate->setValue(testCaseSpec.r);
1400	qRate->setValue(testCaseSpec.q);
1401	vol->setValue(testCaseSpec.volatility);
1402
1403	VanillaOption option(
1404	ext::make_shared<PlainVanillaPayoff>(
1405	args: testCaseSpec.optionType, args: testCaseSpec.strike),
1406	ext::make_shared<AmericanExercise>(args: today, args: maturityDate)
1407	);
1408	option.setPricingEngine(qrPlusAmericanEngine);
1409
1410	const Real calculated = option.NPV();
1411	const Real expected = testCaseSpec.expectedValue;
1412
1413	if ((testCaseSpec.precision > 0
1414	&& std::abs(x: expected-calculated) > testCaseSpec.precision)
1415	|| (testCaseSpec.precision < 0
1416	&& expected > 0.1 && std::abs(x: calculated-expected)/expected > 0.005)
1417	|| (testCaseSpec.precision < 0 && expected <= 0.1
1418	&& std::abs(x: expected-calculated) > 5e-4)) {
1419	BOOST_ERROR("QR+ boundary approximation failed to "
1420	"reproduce cached edge and PDE values for "
1421	<< "\n OptionType: " <<
1422	((testCaseSpec.optionType == Option::Call)? "Call" : "Put")
1423	<< std::setprecision(16)
1424	<< "\n spot: " << spot->value()
1425	<< "\n strike: " << testCaseSpec.strike
1426	<< "\n r: " << rRate->value()
1427	<< "\n q: " << qRate->value()
1428	<< "\n vol: " << vol->value()
1429	<< "\n calculated: " << calculated
1430	<< "\n expected: " << expected);
1431	}
1432	};
1433	}
1434
1435	void AmericanOptionTest::testQdFpIterationScheme() {
1436	BOOST_TEST_MESSAGE("Testing Legendre and tanh-sinh iteration "
1437	"scheme for QD+ fixed-point American engine...");
1438
1439	const Real tol = 1e-8;
1440	const Size l=32, m=6, n=18, p=36;
1441
1442	const ext::shared_ptr<QdFpIterationScheme> schemes[] = {
1443	ext::make_shared<QdFpLegendreScheme>(args: l, args: m, args: n, args: p),
1444	ext::make_shared<QdFpLegendreTanhSinhScheme>(args: l, args: m, args: n, args: tol),
1445	ext::make_shared<QdFpTanhSinhIterationScheme>(args: m, args: n, args: tol)
1446	};
1447
1448	const NormalDistribution nd;
1449
1450	for (const auto& scheme: schemes) {
1451	BOOST_CHECK_EQUAL(n, scheme->getNumberOfChebyshevInterpolationNodes());
1452	BOOST_CHECK_EQUAL(1, scheme->getNumberOfJacobiNewtonFixedPointSteps());
1453	BOOST_CHECK_EQUAL(m-1, scheme->getNumberOfNaiveFixedPointSteps());
1454
1455	QL_CHECK_SMALL(scheme->getFixedPointIntegrator()
1456	->operator()(nd, -10.0, 10.0) - 1.0, tol);
1457	QL_CHECK_SMALL(scheme->getExerciseBoundaryToPriceIntegrator()
1458	->operator()(nd, -10.0, 10.0) - 1.0, tol);
1459	}
1460	}
1461
1462
1463	void AmericanOptionTest::testAndersenLakeHighPrecisionExample() {
1464	BOOST_TEST_MESSAGE("Testing Andersen, Lake and Offengenden "
1465	"high precision example...");
1466
1467	// Example and results are taken from
1468	// Leif Andersen, Mark Lake and Dimitri Offengenden (2015)
1469	// "High Performance American Option Pricing",
1470	// https://papers.ssrn.com/sol3/papers.cfm?abstract_id=2547027
1471
1472	struct SchemeSpec {
1473	Size l, m, n;
1474	Rate r;
1475	Real expected[2];
1476	Real tol;
1477	};
1478
1479	const SchemeSpec testCases[] = {
1480	{ .l: 24, .m: 3, .n: 9, .r: 0.05, .expected: {0.1069528125898476, 0.1069524359360852}, .tol: 1e-6},
1481	{ .l: 5, .m: 1, .n: 4, .r: 0.05, .expected: {0.1070237787625299, 0.1070042740171235}, .tol: 1e-3},
1482	{ .l: 11, .m: 2, .n: 5, .r: 0.05, .expected: {0.106938750864602, 0.1069479057531648}, .tol: 1e-4},
1483	{ .l: 35, .m: 8, .n: 16, .r: 0.05, .expected: {0.1069527032381714, 0.106952558361499}, .tol: 1e-9},
1484	{ .l: 65, .m: 8, .n: 32, .r: 0.05, .expected: {0.1069527028247546, 0.1069526779971959}, .tol: 1e-11},
1485	{ .l: 5, .m: 1, .n: 4, .r: 0.075, .expected: {0.3674420299196104, 0.3674766444325588}, .tol: 1e-3},
1486	{ .l: 11, .m: 2, .n: 5, .r: 0.075, .expected: {0.3671056766787473, 0.3671024005532715}, .tol: 1e-4},
1487	{ .l: 35, .m: 8, .n: 16,.r: 0.075, .expected: {0.3671116758420414, 0.3671111055677869}, .tol: 1e-9},
1488	{ .l: 65, .m: 8, .n: 32,.r: 0.075, .expected: {0.3671112309062572, 0.3671111267813689}, .tol: 1e-11}
1489	};
1490
1491	const DayCounter dc = Actual365Fixed();
1492	const Date today = Date(25, July, 2022);
1493	Settings::instance().evaluationDate() = today;
1494
1495	const auto spot = ext::make_shared<SimpleQuote>(args: 100.0);
1496	const Real strike = 100.0;
1497	const Rate q = 0.05;
1498	const Volatility vol = 0.25;
1499	const Date maturityDate = today + Period(1, Years);
1500
1501	const auto payoff =
1502	ext::make_shared<PlainVanillaPayoff>(args: Option::Put, args: strike);
1503
1504	for (const auto& testCase: testCases) {
1505	const Size l = testCase.l;
1506	const Size m = testCase.m;
1507	const Size n = testCase.n;
1508	const Rate r = testCase.r;
1509	const Real tol = testCase.tol;
1510
1511	const auto bsProcess = ext::make_shared<BlackScholesMertonProcess>(
1512	args: Handle<Quote>(spot),
1513	args: Handle<YieldTermStructure>(flatRate(today, forward: q, dc)),
1514	args: Handle<YieldTermStructure>(flatRate(today, forward: r, dc)),
1515	args: Handle<BlackVolTermStructure>(flatVol(today, volatility: vol, dc))
1516	);
1517
1518	VanillaOption americanOption(
1519	payoff, ext::make_shared<AmericanExercise>(args: today, args: maturityDate));
1520
1521	VanillaOption europeanOption(
1522	payoff, ext::make_shared<EuropeanExercise>(args: maturityDate));
1523
1524	europeanOption.setPricingEngine(
1525	ext::make_shared<AnalyticEuropeanEngine>(args: bsProcess));
1526
1527	const Real europeanNPV = europeanOption.NPV();
1528
1529	const QdFpAmericanEngine::FixedPointEquation schemes[] = {
1530	QdFpAmericanEngine::FP_A, QdFpAmericanEngine::FP_B
1531	};
1532
1533	for (Size i=0; i < LENGTH(schemes); ++i) {
1534
1535	americanOption.setPricingEngine(
1536	ext::make_shared<QdFpAmericanEngine>(
1537	args: bsProcess,
1538	args: ext::make_shared<QdFpLegendreTanhSinhScheme>(args: l, args: m, args: n, args: tol),
1539	args: schemes[i])
1540	);
1541
1542	const Real americanNPV = americanOption.NPV();
1543	const Real americanPremium = americanNPV - europeanNPV;
1544
1545	const Real diff = std::abs(x: americanPremium - testCase.expected[i]);
1546	if (diff > tol) {
1547	BOOST_ERROR("failed to reproduce high precision literature values"
1548	<< "\n FP-Scheme: " <<
1549	((schemes[i] == QdFpAmericanEngine::FP_A)? "FP-A" : "FP-B")
1550	<< "\n r : " << r
1551	<< "\n (l,m,n) : (" << l << "," << m << "," << n << ")"
1552	<< "\n diff : " << diff
1553	<< "\n tol : " << tol);
1554	}
1555	}
1556	}
1557	}
1558
1559
1560	void AmericanOptionTest::testQdEngineStandardExample() {
1561	BOOST_TEST_MESSAGE("Testing Andersen, Lake and Offengenden "
1562	"standard example...");
1563
1564	const DayCounter dc = Actual365Fixed();
1565	const Date today = Date(1, June, 2022);
1566	Settings::instance().evaluationDate() = today;
1567
1568	const Real S = 100;
1569	const Real K = 95;
1570	const Rate r = 0.075;
1571	const Rate q = 0.05;
1572	const Volatility sigma = 0.25;
1573	const Date maturityDate = today + Period(1, Years);
1574
1575	const auto bsProcess = ext::make_shared<BlackScholesMertonProcess>(
1576	args: Handle<Quote>(ext::make_shared<SimpleQuote>(args: S)),
1577	args: Handle<YieldTermStructure>(flatRate(today, forward: q, dc)),
1578	args: Handle<YieldTermStructure>(flatRate(today, forward: r, dc)),
1579	args: Handle<BlackVolTermStructure>(flatVol(today, volatility: sigma, dc))
1580	);
1581
1582	const auto payoff =
1583	ext::make_shared<PlainVanillaPayoff>(args: Option::Put, args: K);
1584
1585	VanillaOption europeanOption(
1586	payoff, ext::make_shared<EuropeanExercise>(args: maturityDate));
1587
1588	europeanOption.setPricingEngine(
1589	ext::make_shared<AnalyticEuropeanEngine>(args: bsProcess));
1590
1591	VanillaOption americanOption(
1592	payoff, ext::make_shared<AmericanExercise>(args: today, args: maturityDate));
1593
1594
1595	const QdFpAmericanEngine::FixedPointEquation schemes[] = {
1596	QdFpAmericanEngine::FP_A, QdFpAmericanEngine::FP_B
1597	};
1598	const Real expected[] = { 0.2386475283369327, 0.2386596962737606 };
1599
1600	for (Size i=0; i < LENGTH(schemes); ++i) {
1601	americanOption.setPricingEngine(
1602	ext::make_shared<QdFpAmericanEngine>(
1603	args: bsProcess,
1604	args: ext::make_shared<QdFpLegendreScheme>(args: 32, args: 2, args: 15, args: 48),
1605	args: schemes[i])
1606	);
1607	const Real calculated = americanOption.NPV() - europeanOption.NPV();
1608
1609	const Real tol = 1e-15;
1610	const Real diff = std::abs(x: calculated - expected[i]);
1611
1612	if (diff > tol) {
1613	BOOST_ERROR("failed to reproduce high precision test values"
1614	<< "\n diff : " << diff
1615	<< "\n tol : " << tol);
1616	}
1617	}
1618	}
1619
1620	namespace {
1621	class QdFpGaussLobattoScheme: public QdFpIterationScheme {
1622	public:
1623	QdFpGaussLobattoScheme(Size m, Size n, Real eps)
1624	: m_(m), n_(n),
1625	integrator_(ext::make_shared<GaussLobattoIntegral>(
1626	args: 100000, QL_MAX_REAL, args: 0.1*eps)) {
1627	}
1628	Size getNumberOfChebyshevInterpolationNodes() const override {
1629	return n_;
1630	}
1631	Size getNumberOfNaiveFixedPointSteps() const override {
1632	return m_-1;
1633	}
1634	Size getNumberOfJacobiNewtonFixedPointSteps() const override {
1635	return Size(1);
1636	}
1637	ext::shared_ptr<Integrator>
1638	getFixedPointIntegrator() const override {
1639	return integrator_;
1640	}
1641	ext::shared_ptr<Integrator>
1642	getExerciseBoundaryToPriceIntegrator() const override {
1643	return integrator_;
1644	}
1645
1646	private:
1647	const Size m_, n_;
1648	const ext::shared_ptr<Integrator> integrator_;
1649	};
1650	}
1651
1652	void AmericanOptionTest::testBulkQdFpAmericanEngine() {
1653	BOOST_TEST_MESSAGE("Testing Andersen, Lake and Offengenden "
1654	"bulk examples...");
1655
1656	// Examples are taken from
1657	// Leif Andersen, Mark Lake and Dimitri Offengenden (2015)
1658	// "High Performance American Option Pricing",
1659	// https://papers.ssrn.com/sol3/papers.cfm?abstract_id=2547027
1660
1661	const DayCounter dc = Actual365Fixed();
1662	const Date today = Date(1, June, 2022);
1663	Settings::instance().evaluationDate() = today;
1664
1665	const auto spot = ext::make_shared<SimpleQuote>(args: 1.0);
1666	const auto rRate = ext::make_shared<SimpleQuote>(args: 0.0);
1667	const auto qRate = ext::make_shared<SimpleQuote>(args: 0.0);
1668	const auto vol = ext::make_shared<SimpleQuote>(args: 0.0);
1669
1670	// original test set from the article, takes too long
1671	// const Size T[] = {30, 91, 182, 273, 365};
1672	// const Rate rf[] = {0.02, 0.04, 0.06, 0.08, 0.1};
1673	// const Rate qy[] = {0, 0.04, 0.08, 0.12};
1674	// const Real S[] = {25, 50, 80, 90, 100, 110, 120, 150, 175, 200};
1675	// const Volatility sig[] = {0.1, 0.2, 0.3, 0.4, 0.5, 0.6};
1676
1677	const Size T[] = {30, 182, 365};
1678	const Rate rf[] = {0.02, 0.04, 0.06, 0.1};
1679	const Rate qy[] = {0, 0.04, 0.08, 0.12};
1680	const Real S[] = {25, 75, 100, 125, 200};
1681	const Volatility sig[] = {0.1, 0.25, 0.6};
1682
1683	const auto payoff = ext::make_shared<PlainVanillaPayoff>(args: Option::Put, args: 100);
1684
1685	const auto bsProcess = ext::make_shared<BlackScholesMertonProcess>(
1686	args: Handle<Quote>(spot),
1687	args: Handle<YieldTermStructure>(flatRate(today, forward: qRate, dc)),
1688	args: Handle<YieldTermStructure>(flatRate(today, forward: rRate, dc)),
1689	args: Handle<BlackVolTermStructure>(flatVol(today, volatility: vol, dc))
1690	);
1691
1692	const auto qdFpFastAmericanEngine =
1693	ext::make_shared<QdFpAmericanEngine>(
1694	args: bsProcess, args: QdFpAmericanEngine::fastScheme());
1695
1696	const auto qdFpAccurateAmericanEngine =
1697	ext::make_shared<QdFpAmericanEngine>(
1698	args: bsProcess, args: QdFpAmericanEngine::accurateScheme());
1699
1700	const auto qdFpGaussLobattoAmericanEngine =
1701	ext::make_shared<QdFpAmericanEngine>(
1702	args: bsProcess,args: ext::make_shared<QdFpGaussLobattoScheme>(args: 3, args: 7, args: 1e-5));
1703
1704	IncrementalStatistics statsAccurate, statsLobatto;
1705	for (auto t: T) {
1706	const Date maturityDate = today + Period(t, Days);
1707	VanillaOption option(
1708	payoff, ext::make_shared<AmericanExercise>(args: today, args: maturityDate));
1709
1710	for (auto r: rf) {
1711	rRate->setValue(r);
1712	for (auto q: qy) {
1713	qRate->setValue(q);
1714	for (auto v: sig) {
1715	vol->setValue(v);
1716	for (auto s: S) {
1717	spot->setValue(s);
1718
1719	option.setPricingEngine(qdFpFastAmericanEngine);
1720	const Real fast = option.NPV();
1721
1722	option.setPricingEngine(qdFpAccurateAmericanEngine);
1723	const Real accurate = option.NPV();
1724	statsAccurate.add(value: std::abs(x: fast-accurate));
1725
1726	option.setPricingEngine(qdFpGaussLobattoAmericanEngine);
1727	const Real lobatto = option.NPV();
1728	statsLobatto.add(value: std::abs(x: accurate-lobatto));
1729	}
1730	}
1731	}
1732	}
1733	}
1734
1735
1736	const Real tolStdDev = 1e-4;
1737	if (statsAccurate.standardDeviation() > tolStdDev)
1738	BOOST_ERROR("failed to reproduce low RMSE with fast American engine"
1739	<< "\n RMSE diff: " << statsAccurate.standardDeviation()
1740	<< "\n tol : " << tolStdDev);
1741
1742	if (statsLobatto.standardDeviation() > tolStdDev)
1743	BOOST_ERROR("failed to reproduce low RMSE with fast Lobatto American engine"
1744	<< "\n RMSE diff: " << statsLobatto.standardDeviation()
1745	<< "\n tol : " << tolStdDev);
1746
1747	const Real tolMax = 2.5e-3;
1748	if (statsAccurate.max() > tolMax)
1749	BOOST_ERROR("failed to reproduce low max deviation "
1750	"with fast American engine"
1751	<< "\n max diff: " << statsAccurate.max()
1752	<< "\n tol : " << tolMax);
1753
1754	if (statsLobatto.max() > tolMax)
1755	BOOST_ERROR("failed to reproduce low max deviation "
1756	"with fast Lobatto American engine"
1757	<< "\n max diff: " << statsLobatto.max()
1758	<< "\n tol : " << tolMax);
1759	}
1760
1761	void AmericanOptionTest::testQdEngineWithLobattoIntegral() {
1762	BOOST_TEST_MESSAGE("Testing Andersen, Lake and Offengenden "
1763	"with high precision Gauss-Lobatto integration...");
1764
1765	const DayCounter dc = Actual365Fixed();
1766	const Date today = Date(5, November, 2022);
1767	Settings::instance().evaluationDate() = today;
1768
1769	const auto spot = ext::make_shared<SimpleQuote>(args: 36);
1770	const Real K = 40;
1771	const Rate r = 0.075;
1772	const Rate q = 0.01;
1773	const Volatility sigma = 0.40;
1774	const Date maturityDate = today + Period(2, Years);
1775
1776	const auto bsProcess = ext::make_shared<BlackScholesMertonProcess>(
1777	args: Handle<Quote>(spot),
1778	args: Handle<YieldTermStructure>(flatRate(today, forward: q, dc)),
1779	args: Handle<YieldTermStructure>(flatRate(today, forward: r, dc)),
1780	args: Handle<BlackVolTermStructure>(flatVol(today, volatility: sigma, dc))
1781	);
1782
1783	VanillaOption option(
1784	ext::make_shared<PlainVanillaPayoff>(args: Option::Put, args: K),
1785	ext::make_shared<AmericanExercise>(args: today, args: maturityDate)
1786	);
1787
1788	const QdFpAmericanEngine::FixedPointEquation schemes[] = {
1789	QdFpAmericanEngine::FP_A, QdFpAmericanEngine::FP_B
1790	};
1791
1792	const auto gaussLobattoScheme =
1793	ext::make_shared<QdFpGaussLobattoScheme>(args: 10, args: 30, args: 1e-10);
1794
1795	for (auto scheme: schemes) {
1796	const auto highPrecisionEngine =
1797	ext::make_shared<QdFpAmericanEngine>(
1798	args: bsProcess, args: QdFpAmericanEngine::highPrecisionScheme(), args&: scheme);
1799	const auto lobattoEngine =
1800	ext::make_shared<QdFpAmericanEngine>(
1801	args: bsProcess, args: gaussLobattoScheme, args&: scheme);
1802
1803	for (Real s: std::list<Real>{36, 40-1e-8, 40, 40+1e-8, 50}) {
1804	spot->setValue(s);
1805
1806	option.setPricingEngine(highPrecisionEngine);
1807	const Real highPrecisionNPV = option.NPV();
1808
1809	option.setPricingEngine(lobattoEngine);
1810	const Real lobattoNPV = option.NPV();
1811
1812	const Real diff = std::abs(x: lobattoNPV - highPrecisionNPV);
1813	const Real tol = 1e-11;
1814
1815	if (diff > tol || std::isnan(x: lobattoNPV)) {
1816	BOOST_ERROR("failed to reproduce high precision American "
1817	"option values with QD+ fixed point and Lobatto integration"
1818	<< "\n FP-Scheme: " <<
1819	((scheme == QdFpAmericanEngine::FP_A)? "FP-A" : "FP-B")
1820	<< "\n spot : " << s
1821	<< "\n diff : " << diff
1822	<< "\n tol : " << tol);
1823	}
1824	}
1825	}
1826	}
1827
1828	void AmericanOptionTest::testQdNegativeDividendYield() {
1829	BOOST_TEST_MESSAGE("Testing Andersen, Lake and Offengenden "
1830	"with positive or zero interest rate and "
1831	"negative dividend yield...");
1832
1833	const DayCounter dc = Actual365Fixed();
1834	const Date today = Date(5, December, 2022);
1835	Settings::instance().evaluationDate() = today;
1836	const Date maturityDate = today + Period(18, Months);
1837
1838	const Real K = 90;
1839	const auto spot = ext::make_shared<SimpleQuote>(args: 100);
1840	const auto qRate = ext::make_shared<SimpleQuote>(args: 0);
1841	const auto rRate = ext::make_shared<SimpleQuote>(args: 0);
1842
1843	const auto process = ext::make_shared<BlackScholesMertonProcess>(
1844	args: Handle<Quote>(spot),
1845	args: Handle<YieldTermStructure>(flatRate(forward: qRate, dc)),
1846	args: Handle<YieldTermStructure>(flatRate(forward: rRate, dc)),
1847	args: Handle<BlackVolTermStructure>(flatVol(volatility: 0.4, dc))
1848	);
1849
1850	VanillaOption option(
1851	ext::make_shared<PlainVanillaPayoff>(args: Option::Put, args: K),
1852	ext::make_shared<AmericanExercise>(args: today, args: maturityDate)
1853	);
1854	const auto qdPlusEngine =
1855	ext::make_shared<QdPlusAmericanEngine>(args: process);
1856	const auto qdFpEngine =
1857	ext::make_shared<QdFpAmericanEngine>(args: process);
1858	const auto fdmEngine =
1859	ext::make_shared<FdBlackScholesVanillaEngine>(args: process, args: 800, args: 200);
1860
1861	const Rate rRates[] = {0.025, 0.5, 0.0, 1e-6};
1862	const Rate qRates[] = {-0.05, -0.1, -0.5, 0.0};
1863
1864	for (auto r: rRates) {
1865	rRate->setValue(r);
1866	for (auto q: qRates) {
1867	qRate->setValue(q);
1868
1869	option.setPricingEngine(qdFpEngine);
1870	const Real qdFpNPV = option.NPV();
1871	option.setPricingEngine(qdPlusEngine);
1872	const Real qdPlusNPV = option.NPV();
1873	option.setPricingEngine(fdmEngine);
1874	const Real fdmNPV = option.NPV();
1875
1876	const Real tol = 1.5e-2;
1877	const Real diffFdmQqFp = std::abs(x: fdmNPV - qdFpNPV);
1878
1879	if (diffFdmQqFp > tol || std::isnan(x: diffFdmQqFp)) {
1880	BOOST_ERROR("failed to reproduce QD+ fixed point values"
1881	<< "\n r : " << r
1882	<< "\n q : " << q
1883	<< "\n diff : " << diffFdmQqFp
1884	<< "\n tol : " << tol);
1885	}
1886
1887	const Real diffFdmQdPlus = std::abs(x: fdmNPV - qdPlusNPV);
1888	if (diffFdmQdPlus > 5*tol || std::isnan(x: diffFdmQdPlus)) {
1889	BOOST_ERROR("failed to reproduce QD+ values"
1890	<< "\n r : " << r
1891	<< "\n q : " << q
1892	<< "\n diff : " << diffFdmQdPlus
1893	<< "\n tol : " << 5*tol);
1894	}
1895	}
1896	}
1897	}
1898
1899	void AmericanOptionTest::testBjerksundStenslandEuropeanGreeks() {
1900	BOOST_TEST_MESSAGE("Testing Bjerksund-Stensland greeks when early "
1901	"exercise is not optimal...");
1902
1903	const Date today = Date(5, November, 2022);
1904	Settings::instance().evaluationDate() = today;
1905
1906	const auto spot = ext::make_shared<SimpleQuote>(args: 100);
1907	const Real K = 105;
1908
1909	const Volatility sigma = 0.40;
1910	const Date maturityDate = today + Period(724, Days);
1911
1912	const auto qTS = ext::make_shared<SimpleQuote>(args: 0.0);
1913	const auto rTS = ext::make_shared<SimpleQuote>(args: 0.0);
1914
1915	const auto bsProcess = ext::make_shared<BlackScholesMertonProcess>(
1916	args: Handle<Quote>(spot),
1917	args: Handle<YieldTermStructure>(flatRate(forward: qTS, dc: Actual365Fixed())),
1918	args: Handle<YieldTermStructure>(flatRate(forward: rTS, dc: Actual360())),
1919	args: Handle<BlackVolTermStructure>(
1920	flatVol(today, volatility: sigma, dc: Thirty360(Thirty360::European)))
1921	);
1922
1923	struct OptionSpec {
1924	Option::Type type;
1925	Real r;
1926	Real q;
1927	};
1928
1929	const OptionSpec testCaseSpecs[] = {
1930	{.type: Option::Put, .r: -0.05, .q: 0.02},
1931	{.type: Option::Call, .r: 0.05, .q: -0.025}
1932	};
1933
1934	const auto europeanExercise =
1935	ext::make_shared<EuropeanExercise>(args: maturityDate);
1936	const auto americanExercise =
1937	ext::make_shared<AmericanExercise>(args: today, args: maturityDate);
1938
1939	const auto europeanEngine =
1940	ext::make_shared<AnalyticEuropeanEngine>(args: bsProcess);
1941
1942	const auto bjerksundStenslandEngine =
1943	ext::make_shared<BjerksundStenslandApproximationEngine>(args: bsProcess);
1944
1945
1946	for (const auto& testCaseSpec: testCaseSpecs) {
1947	qTS->setValue(testCaseSpec.q);
1948	rTS->setValue(testCaseSpec.r);
1949
1950	VanillaOption americanOption(
1951	ext::make_shared<PlainVanillaPayoff>(args: testCaseSpec.type, args: K),
1952	americanExercise
1953	);
1954	americanOption.setPricingEngine(bjerksundStenslandEngine);
1955
1956	VanillaOption europeanOption(
1957	ext::make_shared<PlainVanillaPayoff>(args: testCaseSpec.type, args: K),
1958	europeanExercise
1959	);
1960	europeanOption.setPricingEngine(europeanEngine);
1961
1962	constexpr double tol = 1000*QL_EPSILON;
1963
1964	BOOST_CHECK_CLOSE(europeanOption.NPV(), americanOption.NPV(), tol);
1965	BOOST_CHECK_CLOSE(europeanOption.delta(), americanOption.delta(), tol);
1966	BOOST_CHECK_CLOSE(europeanOption.strikeSensitivity(), americanOption.strikeSensitivity(), tol);
1967	BOOST_CHECK_CLOSE(europeanOption.gamma(), americanOption.gamma(), tol);
1968	BOOST_CHECK_CLOSE(europeanOption.vega(), americanOption.vega(), tol);
1969	BOOST_CHECK_CLOSE(europeanOption.theta(), americanOption.theta(), tol);
1970	BOOST_CHECK_CLOSE(europeanOption.thetaPerDay(), americanOption.thetaPerDay(), tol);
1971	BOOST_CHECK_CLOSE(europeanOption.rho(), americanOption.rho(), tol);
1972	BOOST_CHECK_CLOSE(europeanOption.dividendRho(), americanOption.dividendRho(), tol);
1973	}
1974	}
1975
1976
1977	void AmericanOptionTest::testBjerksundStenslandAmericanGreeks() {
1978	BOOST_TEST_MESSAGE("Testing Bjerksund-Stensland American greeks...");
1979
1980	const Date today = Date(5, December, 2022);
1981	Settings::instance().evaluationDate() = today;
1982
1983	const auto spot = ext::make_shared<SimpleQuote>(args: 0);
1984	const auto vol = ext::make_shared<SimpleQuote>(args: 0);
1985
1986	const auto qRate = ext::make_shared<SimpleQuote>(args: 0);
1987	const auto rRate = ext::make_shared<SimpleQuote>(args: 0);
1988
1989	const auto bsProcess = ext::make_shared<BlackScholesMertonProcess>(
1990	args: Handle<Quote>(spot),
1991	args: Handle<YieldTermStructure>(flatRate(forward: qRate, dc: Actual360())),
1992	args: Handle<YieldTermStructure>(flatRate(forward: rRate, dc: Actual365Fixed())),
1993	args: Handle<BlackVolTermStructure>(
1994	flatVol(today, volatility: vol, dc: Thirty360(Thirty360::ISDA)))
1995	);
1996
1997	const auto bjerksundStenslandEngine =
1998	ext::make_shared<BjerksundStenslandApproximationEngine>(args: bsProcess);
1999
2000	const Real strike = 100;
2001	const Option::Type types[] = { Option::Call, Option::Put };
2002	const Rate rf[] = {0.0, 0.02, 0.06, 0.1, 0.2};
2003	const Rate qy[] = {0.0, 0.08, 0.12};
2004
2005	const Volatility sig[] = {0.1, 0.2, 0.4, 1.0};
2006	const Real S[] = {25, 50, 99.9, 110, 150, 200};
2007	const Size T[] = {30, 182, 365, 1825};
2008
2009	const Real f_d = 1e-5;
2010	const Real f_g = 5e-5;
2011	const Real f_q = 1e-6;
2012
2013	for (auto type: types) {
2014	const auto payoff = [type](Real strike) {
2015	return ext::make_shared<PlainVanillaPayoff>(args: type, args&: strike);};
2016
2017	const auto stdPayoff = payoff(strike);
2018
2019	for (auto t: T) {
2020	const Date maturityDate = today + Period(t, Days);
2021	const auto exercise = [today, maturityDate](const Period& offset) {
2022	return ext::make_shared<AmericanExercise>(
2023	args: today, args: maturityDate + offset);};
2024
2025	const auto stdExercise = exercise(Period(0, Days));
2026
2027	VanillaOption option(stdPayoff, stdExercise);
2028	option.setPricingEngine(bjerksundStenslandEngine);
2029
2030	VanillaOption strike_up(payoff(strike*(1+f_d)), stdExercise);
2031	strike_up.setPricingEngine(bjerksundStenslandEngine);
2032	VanillaOption strike_down(payoff(strike*(1-f_d)), stdExercise);
2033	strike_down.setPricingEngine(bjerksundStenslandEngine);
2034
2035	VanillaOption day_up(stdPayoff, exercise(Period(1, Days)));
2036	day_up.setPricingEngine(bjerksundStenslandEngine);
2037	VanillaOption day_down(stdPayoff, exercise(Period(-1, Days)));
2038	day_down.setPricingEngine(bjerksundStenslandEngine);
2039
2040
2041	for (auto r: rf) {
2042	rRate->setValue(r);
2043	for (auto q: qy) {
2044	qRate->setValue(q);
2045	for (auto v: sig) {
2046	vol->setValue(v);
2047	for (auto s: S) {
2048	spot->setValue(s);
2049
2050	const Real npv = option.NPV();
2051	const Real delta = option.delta();
2052	const Real gamma = option.gamma();
2053	const Real strikeSensitivity = option.strikeSensitivity();
2054	const Real dividendRho = option.dividendRho();
2055	const Real rho = option.rho();
2056	const Real vega = option.vega();
2057	const Real theta = option.theta();
2058	const std::string exerciseType = ext::any_cast<std::string>(
2059	operand: option.additionalResults().find(x: "exerciseType")->second);
2060
2061	OneAssetOption::results numericalResults;
2062
2063	spot->setValue(s*(1+f_d));
2064	const Real f2 = option.NPV();
2065	spot->setValue(s*(1-f_d));
2066	const Real f1 = option.NPV();
2067	spot->setValue(s);
2068	numericalResults.delta = (f2 - f1)/(2*f_d*s);
2069
2070	Real error = std::abs(x: delta - numericalResults.delta);
2071	if (error > 5e-6)
2072	REPORT_FAILURE("delta", \
2073	stdPayoff, stdExercise, s, q, r, today, v, \
2074	numericalResults.delta, delta, error, 5e-6);
2075
2076	spot->setValue(s*(1+2*f_g));
2077	const Real gp2 = option.NPV();
2078	spot->setValue(s*(1+f_g));
2079	const Real gp1 = option.NPV();
2080	spot->setValue(s*(1-f_g));
2081	const Real gm1 = option.NPV();
2082	spot->setValue(s*(1-2*f_g));
2083	const Real gm2 = option.NPV();
2084	spot->setValue(s);
2085	numericalResults.gamma
2086	= (-gp2 + 16*gp1 - 30*npv + 16*gm1 - gm2)/(12*squared(x: f_g*s));
2087
2088	error = std::abs(x: gamma - numericalResults.gamma);
2089	if (error > 1e-4 && t < 1000)
2090	REPORT_FAILURE("gamma", \
2091	stdPayoff, stdExercise, s, q, r, today, v, \
2092	numericalResults.gamma, gamma, error, 5e-5);
2093
2094	const Real k2 = strike_up.NPV();
2095	const Real k1 = strike_down.NPV();
2096	numericalResults.strikeSensitivity = (k2 - k1)/(2*f_d*strike);
2097	error = std::abs(x: strikeSensitivity - numericalResults.strikeSensitivity);
2098
2099	if (error > 5e-6)
2100	REPORT_FAILURE("strikeSensitivity", \
2101	stdPayoff, stdExercise, s, q, r, today, v, \
2102	numericalResults.strikeSensitivity, strikeSensitivity, error, 5e-6);
2103
2104	if (q != 0.0) {
2105	qRate->setValue(q + f_q);
2106	const Real q2 = option.NPV();
2107	qRate->setValue(q - f_q);
2108	const Real q1 = option.NPV();
2109	qRate->setValue(q);
2110	numericalResults.dividendRho = (q2-q1)/(2*f_q);
2111
2112	error = std::abs(x: dividendRho - numericalResults.dividendRho);
2113
2114	if (error > 3e-2)
2115	REPORT_FAILURE("dividendRho", \
2116	stdPayoff, stdExercise, s, q, r, today, v, \
2117	numericalResults.dividendRho, dividendRho, error, 1e-3);
2118
2119	rRate->setValue(r + f_q);
2120	const Real r2 = option.NPV();
2121	rRate->setValue(r - f_q);
2122	const Real r1 = option.NPV();
2123	rRate->setValue(r);
2124	numericalResults.rho = (r2 - r1)/(2*f_q);
2125
2126	error = std::abs(x: rho - numericalResults.rho);
2127	if (error > 3e-2)
2128	REPORT_FAILURE("rho", \
2129	stdPayoff, stdExercise, s, q, r, today, v, \
2130	numericalResults.rho, rho, error, 1e-3);
2131	}
2132
2133	vol->setValue(v + f_d);
2134	const Real v2 = option.NPV();
2135	vol->setValue(v - f_d);
2136	const Real v1 = option.NPV();
2137	vol->setValue(v);
2138	numericalResults.vega = (v2 - v1)/(2*f_d);
2139
2140	error = std::abs(x: vega - numericalResults.vega);
2141	if (error > 5e-4)
2142	REPORT_FAILURE("vega", \
2143	stdPayoff, stdExercise, s, q, r, today, v, \
2144	numericalResults.vega, vega, error, 5e-4);
2145
2146	if (exerciseType == "American") {
2147	const Real t2 = day_up.NPV();
2148	const Real t1 = day_down.NPV();
2149	numericalResults.thetaPerDay = (t1-t2)/2;
2150	numericalResults.theta = 365*numericalResults.thetaPerDay;
2151	error = std::abs(x: theta - numericalResults.theta);
2152	const Real thetaTol = (t < 60) ? 3.0: 5e-4;
2153	if (error > thetaTol) {
2154	REPORT_FAILURE("theta", \
2155	stdPayoff, stdExercise, s, q, r, today, v, \
2156	numericalResults.theta, theta, error, thetaTol);
2157	}
2158	}
2159	}
2160	}
2161	}
2162	}
2163	}
2164	}
2165	}
2166
2167
2168	void AmericanOptionTest::testSingleBjerksundStenslandGreeks() {
2169	BOOST_TEST_MESSAGE("Testing a single Bjerksund-Stensland greeks set...");
2170
2171	const Date today = Date(20, January, 2023);
2172	Settings::instance().evaluationDate() = today;
2173
2174	const Real s = 100;
2175	const Volatility v = 0.3;
2176	const Rate q = 0.04;
2177	const Rate r = 0.07;
2178
2179	const auto spot = ext::make_shared<SimpleQuote>(args: s);
2180	const auto vol = ext::make_shared<SimpleQuote>(args: v);
2181
2182	const auto qRate = ext::make_shared<SimpleQuote>(args: q);
2183	const auto rRate = ext::make_shared<SimpleQuote>(args: r);
2184
2185	const auto bsProcess = ext::make_shared<BlackScholesMertonProcess>(
2186	args: Handle<Quote>(spot),
2187	args: Handle<YieldTermStructure>(flatRate(forward: qRate, dc: Actual365Fixed())),
2188	args: Handle<YieldTermStructure>(flatRate(forward: rRate, dc: Actual365Fixed())),
2189	args: Handle<BlackVolTermStructure>(
2190	flatVol(today, volatility: vol, dc: Actual365Fixed()))
2191	);
2192
2193	const Date maturityDate = today + Period(2, Years);
2194
2195	const auto exercise
2196	= ext::make_shared<AmericanExercise>(args: today, args: maturityDate);
2197	const auto payoff
2198	= ext::make_shared<PlainVanillaPayoff>(args: Option::Call, args: 100);
2199
2200	VanillaOption option(payoff, exercise);
2201
2202	option.setPricingEngine(
2203	ext::make_shared<BjerksundStenslandApproximationEngine>(args: bsProcess)
2204	);
2205
2206	const Real npv = option.NPV();
2207	const Real delta = option.delta();
2208	const Real gamma = option.gamma();
2209	const Real strikeSensitivity = option.strikeSensitivity();
2210	const Real divRho = option.dividendRho();
2211	const Real rho = option.rho();
2212	const Real vega = option.vega();
2213	const Real theta = option.theta();
2214	const Real thetaPerDay = option.thetaPerDay();
2215	const std::string exerciseType = ext::any_cast<std::string>(
2216	operand: option.additionalResults().find(x: "exerciseType")->second);
2217
2218	const Real expectedNpv = 17.9251834488399169;
2219	const Real expectedDelta = 0.590801845261082592;
2220	const Real expectedGamma = 0.00825347110063545664;
2221	const Real expectedStrikeSensitivity = -0.411550010772683383;
2222	const Real expectedDivRho = -114.137818682236826;
2223	const Real expectedRho = 80.4900013901554416;
2224	const Real expectedVega = 49.2906331545933227;
2225	const Real expectedTheta = -4.22540293840206704;
2226
2227	const auto report = [=](Real value, Real expectedValue, const std::string& name) {
2228	constexpr double tol = 1e6*QL_EPSILON;
2229	const Real error = std::abs(x: value-expectedValue);
2230	if (error > tol)
2231	REPORT_FAILURE(name, \
2232	payoff, exercise, s, q, r, today, v, \
2233	value, expectedValue, error, tol);
2234	};
2235
2236	report(npv, expectedNpv, "npv");
2237	report(delta, expectedDelta, "delta");
2238	report(gamma, expectedGamma, "gamma");
2239	report(strikeSensitivity, expectedStrikeSensitivity,
2240	"strikeSensitivity");
2241	report(divRho, expectedDivRho, "dividendRho");
2242	report(rho, expectedRho, "rho");
2243	report(vega, expectedVega, "vega");
2244	report(theta, expectedTheta, "theta");
2245	report(thetaPerDay, expectedTheta/365, "thetaPerDay");
2246
2247	if (exerciseType != "American")
2248	BOOST_FAIL("American exercise type expected");
2249	}
2250
2251
2252	test_suite* AmericanOptionTest::suite(SpeedLevel speed) {
2253	auto* suite = BOOST_TEST_SUITE("American option tests");
2254
2255	suite->add(QUANTLIB_TEST_CASE(&AmericanOptionTest::testBaroneAdesiWhaleyValues));
2256	suite->add(QUANTLIB_TEST_CASE(&AmericanOptionTest::testBjerksundStenslandValues));
2257	suite->add(QUANTLIB_TEST_CASE(&AmericanOptionTest::testJuValues));
2258	suite->add(QUANTLIB_TEST_CASE(&AmericanOptionTest::testFdValues));
2259	suite->add(QUANTLIB_TEST_CASE(&AmericanOptionTest::testFdAmericanGreeks));
2260	suite->add(QUANTLIB_TEST_CASE(&AmericanOptionTest::testFDShoutNPV));
2261	suite->add(QUANTLIB_TEST_CASE(&AmericanOptionTest::testZeroVolFDShoutNPV));
2262	suite->add(QUANTLIB_TEST_CASE(&AmericanOptionTest::testLargeDividendShoutNPV));
2263	suite->add(QUANTLIB_TEST_CASE(&AmericanOptionTest::testEscrowedVsSpotAmericanOption));
2264	suite->add(QUANTLIB_TEST_CASE(&AmericanOptionTest::testTodayIsDividendDate));
2265	suite->add(QUANTLIB_TEST_CASE(&AmericanOptionTest::testCallPutParity));
2266	suite->add(QUANTLIB_TEST_CASE(&AmericanOptionTest::testQdPlusBoundaryValues));
2267	suite->add(QUANTLIB_TEST_CASE(&AmericanOptionTest::testQdPlusBoundaryConvergence));
2268	suite->add(QUANTLIB_TEST_CASE(&AmericanOptionTest::testQdAmericanEngines));
2269	suite->add(QUANTLIB_TEST_CASE(&AmericanOptionTest::testQdFpIterationScheme));
2270	suite->add(QUANTLIB_TEST_CASE(&AmericanOptionTest::testAndersenLakeHighPrecisionExample));
2271	suite->add(QUANTLIB_TEST_CASE(&AmericanOptionTest::testQdEngineStandardExample));
2272	suite->add(QUANTLIB_TEST_CASE(&AmericanOptionTest::testBulkQdFpAmericanEngine));
2273	suite->add(QUANTLIB_TEST_CASE(&AmericanOptionTest::testQdEngineWithLobattoIntegral));
2274	suite->add(QUANTLIB_TEST_CASE(&AmericanOptionTest::testQdNegativeDividendYield));
2275	suite->add(QUANTLIB_TEST_CASE(&AmericanOptionTest::testBjerksundStenslandEuropeanGreeks));
2276	suite->add(QUANTLIB_TEST_CASE(&AmericanOptionTest::testBjerksundStenslandAmericanGreeks));
2277	suite->add(QUANTLIB_TEST_CASE(&AmericanOptionTest::testSingleBjerksundStenslandGreeks));
2278
2279
2280	if (speed <= Fast) {
2281	suite->add(QUANTLIB_TEST_CASE(&AmericanOptionTest::testFdShoutGreeks));
2282	}
2283
2284	return suite;
2285	}
2286
2287