[go: up one dir, main page]
1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2008, 2009, 2014 Klaus Spanderen
5 Copyright (C) 2014 Johannes Göttker-Schnetmann
6
7 This file is part of QuantLib, a free-software/open-source library
8 for financial quantitative analysts and developers - http://quantlib.org/
9
10 QuantLib is free software: you can redistribute it and/or modify it
11 under the terms of the QuantLib license. You should have received a
12 copy of the license along with this program; if not, please email
13 <quantlib-dev@lists.sf.net>. The license is also available online at
14 <http://quantlib.org/license.shtml>.
15
16 This program is distributed in the hope that it will be useful, but WITHOUT
17 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 FOR A PARTICULAR PURPOSE. See the license for more details.
19*/
20
21#include "fdheston.hpp"
22#include "utilities.hpp"
23
24#include <ql/math/functional.hpp>
25#include <ql/quotes/simplequote.hpp>
26#include <ql/time/calendars/target.hpp>
27#include <ql/time/daycounters/actual360.hpp>
28#include <ql/time/daycounters/actualactual.hpp>
29#include <ql/time/daycounters/actual365fixed.hpp>
30#include <ql/instruments/barrieroption.hpp>
31#include <ql/instruments/dividendvanillaoption.hpp>
32#include <ql/models/equity/hestonmodel.hpp>
33#include <ql/termstructures/yield/zerocurve.hpp>
34#include <ql/termstructures/yield/flatforward.hpp>
35#include <ql/termstructures/volatility/equityfx/localconstantvol.hpp>
36#include <ql/methods/finitedifferences/meshers/fdmhestonvariancemesher.hpp>
37#include <ql/pricingengines/barrier/analyticbarrierengine.hpp>
38#include <ql/pricingengines/vanilla/analytichestonengine.hpp>
39#include <ql/pricingengines/vanilla/analyticeuropeanengine.hpp>
40#include <ql/pricingengines/barrier/fdhestonbarrierengine.hpp>
41#include <ql/pricingengines/vanilla/fdhestonvanillaengine.hpp>
42#include <ql/pricingengines/barrier/fdblackscholesbarrierengine.hpp>
43#include <ql/pricingengines/vanilla/fdblackscholesvanillaengine.hpp>
44#include <ql/tuple.hpp>
45
46using namespace QuantLib;
47using boost::unit_test_framework::test_suite;
48
49
50namespace fd_heston_test {
51 struct NewBarrierOptionData {
52 Barrier::Type barrierType;
53 Real barrier;
54 Real rebate;
55 Option::Type type;
56 Real strike;
57 Real s; // spot
58 Rate q; // dividend
59 Rate r; // risk-free rate
60 Time t; // time to maturity
61 Volatility v; // volatility
62 };
63
64 class ParableLocalVolatility : public LocalVolTermStructure {
65 public:
66 ParableLocalVolatility(
67 const Date& referenceDate,
68 Real s0,
69 Real alpha,
70 const DayCounter& dayCounter)
71 : LocalVolTermStructure(
72 referenceDate, NullCalendar(), Following, dayCounter),
73 referenceDate_(referenceDate),
74 s0_(s0),
75 alpha_(alpha) {}
76
77 Date maxDate() const override { return Date::maxDate(); }
78 Real minStrike() const override { return 0.0; }
79 Real maxStrike() const override { return std::numeric_limits<Real>::max(); }
80
81 protected:
82 Volatility localVolImpl(Time, Real s) const override {
83 return alpha_*(squared(x: s0_ - s) + 25.0);
84 }
85
86 private:
87 const Date referenceDate_;
88 const Real s0_, alpha_;
89 };
90}
91
92void FdHestonTest::testFdmHestonVarianceMesher() {
93 BOOST_TEST_MESSAGE("Testing FDM Heston variance mesher...");
94
95 using namespace fd_heston_test;
96
97 const Date today = Date(22, February, 2018);
98 const DayCounter dc = Actual365Fixed();
99 Settings::instance().evaluationDate() = today;
100
101 const ext::shared_ptr<HestonProcess> process(
102 ext::make_shared<HestonProcess>(
103 args: Handle<YieldTermStructure>(flatRate(forward: 0.02, dc)),
104 args: Handle<YieldTermStructure>(flatRate(forward: 0.02, dc)),
105 args: Handle<Quote>(ext::make_shared<SimpleQuote>(args: 100.0)),
106 args: 0.09, args: 1.0, args: 0.09, args: 0.2, args: -0.5));
107
108 const ext::shared_ptr<FdmHestonVarianceMesher> mesher
109 = ext::make_shared<FdmHestonVarianceMesher>(args: 5, args: process, args: 1.0);
110
111 const std::vector<Real> locations = mesher->locations();
112
113 const Real expected[] = {
114 0.0, 6.652314e-02, 9.000000e-02, 1.095781e-01, 2.563610e-01
115 };
116
117 const Real tol = 1e-6;
118 for (Size i=0; i < locations.size(); ++i) {
119 const Real diff = std::fabs(x: expected[i] - locations[i]);
120
121 if (diff > tol) {
122 BOOST_ERROR("Failed to reproduce Heston variance mesh"
123 << "\n calculated: " << locations[i]
124 << "\n expected: " << expected[i]
125 << std::scientific
126 << "\n difference " << diff
127 << "\n tolerance: " << tol);
128 }
129 }
130
131 const ext::shared_ptr<LocalVolTermStructure> lVol =
132 ext::make_shared<LocalConstantVol>(args: today, args: 2.5, args: dc);
133
134 const ext::shared_ptr<FdmHestonLocalVolatilityVarianceMesher> constSlvMesher
135 = ext::make_shared<FdmHestonLocalVolatilityVarianceMesher>
136 (args: 5, args: process, args: lVol, args: 1.0);
137
138 const Real expectedVol = 2.5 * mesher->volaEstimate();
139 const Real calculatedVol = constSlvMesher->volaEstimate();
140
141 const Real diff = std::fabs(x: calculatedVol - expectedVol);
142 if (diff > tol) {
143 BOOST_ERROR("Failed to reproduce Heston local volatility "
144 "variance estimate"
145 << "\n calculated: " << calculatedVol
146 << "\n expected: " << expectedVol
147 << std::scientific
148 << "\n difference " << diff
149 << "\n tolerance: " << tol);
150 }
151
152 const Real alpha = 0.01;
153 const ext::shared_ptr<LocalVolTermStructure> leverageFct
154 = ext::make_shared<ParableLocalVolatility>(args: today, args: 100.0, args: alpha, args: dc);
155
156 const ext::shared_ptr<FdmHestonLocalVolatilityVarianceMesher> slvMesher
157 = ext::make_shared<FdmHestonLocalVolatilityVarianceMesher>(
158 args: 5, args: process, args: leverageFct, args: 0.5, args: 1, args: 0.01);
159
160 const Real initialVolEstimate =
161 ext::make_shared<FdmHestonVarianceMesher>(args: 5, args: process, args: 0.5, args: 1, args: 0.01)->
162 volaEstimate();
163
164 // const Real vEst = leverageFct->localVol(0, 100) * initialVolEstimate;
165 // Mathematica solution
166 // N[Integrate[
167 // alpha*((100*Exp[vEst*x*Sqrt[0.5]] - 100)^2 + 25)*
168 // PDF[NormalDistribution[0, 1], x], {x ,
169 // InverseCDF[NormalDistribution[0, 1], 0.01],
170 // InverseCDF[NormalDistribution[0, 1], 0.99]}]]
171
172 const Real leverageAvg = 0.455881 / (1-0.02);
173
174 const Real volaEstExpected =
175 0.5*(leverageAvg + leverageFct->localVol(t: 0, underlyingLevel: 100)) * initialVolEstimate;
176
177 const Real volaEstCalculated = slvMesher->volaEstimate();
178
179 if (std::fabs(x: volaEstExpected - volaEstCalculated) > 0.001) {
180 BOOST_ERROR("Failed to reproduce Heston local volatility "
181 "variance estimate"
182 << "\n calculated: " << calculatedVol
183 << "\n expected: " << expectedVol
184 << std::scientific
185 << "\n difference " << std::fabs(volaEstExpected - volaEstCalculated)
186 << "\n tolerance: " << tol);
187 }
188}
189
190void FdHestonTest::testFdmHestonBarrierVsBlackScholes() {
191
192 BOOST_TEST_MESSAGE("Testing FDM with barrier option in Heston model...");
193
194 using namespace fd_heston_test;
195
196 NewBarrierOptionData values[] = {
197 /* The data below are from
198 "Option pricing formulas", E.G. Haug, McGraw-Hill 1998 pag. 72
199 */
200 // barrierType, barrier, rebate, type, strike, s, q, r, t, v
201 { .barrierType: Barrier::DownOut, .barrier: 95.0, .rebate: 3.0, .type: Option::Call, .strike: 90, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.25},
202 { .barrierType: Barrier::DownOut, .barrier: 95.0, .rebate: 3.0, .type: Option::Call, .strike: 100, .s: 100.0, .q: 0.00, .r: 0.08, .t: 1.00, .v: 0.30},
203 { .barrierType: Barrier::DownOut, .barrier: 95.0, .rebate: 3.0, .type: Option::Call, .strike: 110, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.25},
204 { .barrierType: Barrier::DownOut, .barrier: 100.0, .rebate: 3.0, .type: Option::Call, .strike: 90, .s: 100.0, .q: 0.00, .r: 0.08, .t: 0.25, .v: 0.25},
205 { .barrierType: Barrier::DownOut, .barrier: 100.0, .rebate: 3.0, .type: Option::Call, .strike: 100, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.25},
206 { .barrierType: Barrier::DownOut, .barrier: 100.0, .rebate: 3.0, .type: Option::Call, .strike: 110, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.25},
207 { .barrierType: Barrier::UpOut, .barrier: 105.0, .rebate: 3.0, .type: Option::Call, .strike: 90, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.25},
208 { .barrierType: Barrier::UpOut, .barrier: 105.0, .rebate: 3.0, .type: Option::Call, .strike: 100, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.25},
209 { .barrierType: Barrier::UpOut, .barrier: 105.0, .rebate: 3.0, .type: Option::Call, .strike: 110, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.25},
210
211 { .barrierType: Barrier::DownIn, .barrier: 95.0, .rebate: 3.0, .type: Option::Call, .strike: 90, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.25},
212 { .barrierType: Barrier::DownIn, .barrier: 95.0, .rebate: 3.0, .type: Option::Call, .strike: 100, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.25},
213 { .barrierType: Barrier::DownIn, .barrier: 95.0, .rebate: 3.0, .type: Option::Call, .strike: 110, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.25},
214 { .barrierType: Barrier::DownIn, .barrier: 100.0, .rebate: 3.0, .type: Option::Call, .strike: 90, .s: 100.0, .q: 0.00, .r: 0.08, .t: 0.25, .v: 0.25},
215 { .barrierType: Barrier::DownIn, .barrier: 100.0, .rebate: 3.0, .type: Option::Call, .strike: 100, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.25},
216 { .barrierType: Barrier::DownIn, .barrier: 100.0, .rebate: 3.0, .type: Option::Call, .strike: 110, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.25},
217 { .barrierType: Barrier::UpIn, .barrier: 105.0, .rebate: 3.0, .type: Option::Call, .strike: 90, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.25},
218 { .barrierType: Barrier::UpIn, .barrier: 105.0, .rebate: 3.0, .type: Option::Call, .strike: 100, .s: 100.0, .q: 0.00, .r: 0.08, .t: 0.40, .v: 0.25},
219 { .barrierType: Barrier::UpIn, .barrier: 105.0, .rebate: 3.0, .type: Option::Call, .strike: 110, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.15},
220
221 { .barrierType: Barrier::DownOut, .barrier: 95.0, .rebate: 3.0, .type: Option::Call, .strike: 90, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.30},
222 { .barrierType: Barrier::DownOut, .barrier: 95.0, .rebate: 3.0, .type: Option::Call, .strike: 100, .s: 100.0, .q: 0.00, .r: 0.08, .t: 0.40, .v: 0.35},
223 { .barrierType: Barrier::DownOut, .barrier: 95.0, .rebate: 3.0, .type: Option::Call, .strike: 110, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.30},
224 { .barrierType: Barrier::DownOut, .barrier: 100.0, .rebate: 3.0, .type: Option::Call, .strike: 90, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.15},
225 { .barrierType: Barrier::DownOut, .barrier: 100.0, .rebate: 3.0, .type: Option::Call, .strike: 100, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.30},
226 { .barrierType: Barrier::DownOut, .barrier: 100.0, .rebate: 3.0, .type: Option::Call, .strike: 110, .s: 100.0, .q: 0.00, .r: 0.00, .t: 1.00, .v: 0.20},
227 { .barrierType: Barrier::UpOut, .barrier: 105.0, .rebate: 3.0, .type: Option::Call, .strike: 90, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.30},
228 { .barrierType: Barrier::UpOut, .barrier: 105.0, .rebate: 3.0, .type: Option::Call, .strike: 100, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.30},
229 { .barrierType: Barrier::UpOut, .barrier: 105.0, .rebate: 3.0, .type: Option::Call, .strike: 110, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.30},
230
231 { .barrierType: Barrier::DownIn, .barrier: 95.0, .rebate: 3.0, .type: Option::Call, .strike: 90, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.30},
232 { .barrierType: Barrier::DownIn, .barrier: 95.0, .rebate: 3.0, .type: Option::Call, .strike: 100, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.30},
233 { .barrierType: Barrier::DownIn, .barrier: 95.0, .rebate: 3.0, .type: Option::Call, .strike: 110, .s: 100.0, .q: 0.00, .r: 0.08, .t: 1.00, .v: 0.30},
234 { .barrierType: Barrier::DownIn, .barrier: 100.0, .rebate: 3.0, .type: Option::Call, .strike: 90, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.30},
235 { .barrierType: Barrier::DownIn, .barrier: 100.0, .rebate: 3.0, .type: Option::Call, .strike: 100, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.30},
236 { .barrierType: Barrier::DownIn, .barrier: 100.0, .rebate: 3.0, .type: Option::Call, .strike: 110, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.30},
237 { .barrierType: Barrier::UpIn, .barrier: 105.0, .rebate: 3.0, .type: Option::Call, .strike: 90, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.30},
238 { .barrierType: Barrier::UpIn, .barrier: 105.0, .rebate: 3.0, .type: Option::Call, .strike: 100, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.30},
239 { .barrierType: Barrier::UpIn, .barrier: 105.0, .rebate: 3.0, .type: Option::Call, .strike: 110, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.30},
240
241 { .barrierType: Barrier::DownOut, .barrier: 95.0, .rebate: 3.0, .type: Option::Put, .strike: 90, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.25},
242 { .barrierType: Barrier::DownOut, .barrier: 95.0, .rebate: 3.0, .type: Option::Put, .strike: 100, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.25},
243 { .barrierType: Barrier::DownOut, .barrier: 95.0, .rebate: 3.0, .type: Option::Put, .strike: 110, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.25},
244 { .barrierType: Barrier::DownOut, .barrier: 100.0, .rebate: 3.0, .type: Option::Put, .strike: 90, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.25},
245 { .barrierType: Barrier::DownOut, .barrier: 100.0, .rebate: 3.0, .type: Option::Put, .strike: 100, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.25},
246 { .barrierType: Barrier::DownOut, .barrier: 100.0, .rebate: 3.0, .type: Option::Put, .strike: 110, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.25},
247 { .barrierType: Barrier::UpOut, .barrier: 105.0, .rebate: 3.0, .type: Option::Put, .strike: 90, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.25},
248 { .barrierType: Barrier::UpOut, .barrier: 105.0, .rebate: 3.0, .type: Option::Put, .strike: 100, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.25},
249 { .barrierType: Barrier::UpOut, .barrier: 105.0, .rebate: 3.0, .type: Option::Put, .strike: 110, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.25},
250
251 { .barrierType: Barrier::DownIn, .barrier: 95.0, .rebate: 3.0, .type: Option::Put, .strike: 90, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.25},
252 { .barrierType: Barrier::DownIn, .barrier: 95.0, .rebate: 3.0, .type: Option::Put, .strike: 100, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.25},
253 { .barrierType: Barrier::DownIn, .barrier: 95.0, .rebate: 3.0, .type: Option::Put, .strike: 110, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.25},
254 { .barrierType: Barrier::DownIn, .barrier: 100.0, .rebate: 3.0, .type: Option::Put, .strike: 90, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.25},
255 { .barrierType: Barrier::DownIn, .barrier: 100.0, .rebate: 3.0, .type: Option::Put, .strike: 100, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.25},
256 { .barrierType: Barrier::DownIn, .barrier: 100.0, .rebate: 3.0, .type: Option::Put, .strike: 110, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.25},
257 { .barrierType: Barrier::UpIn, .barrier: 105.0, .rebate: 3.0, .type: Option::Put, .strike: 90, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.25},
258 { .barrierType: Barrier::UpIn, .barrier: 105.0, .rebate: 3.0, .type: Option::Put, .strike: 100, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.25},
259 { .barrierType: Barrier::UpIn, .barrier: 105.0, .rebate: 3.0, .type: Option::Put, .strike: 110, .s: 100.0, .q: 0.00, .r: 0.04, .t: 1.00, .v: 0.15},
260
261 { .barrierType: Barrier::DownOut, .barrier: 95.0, .rebate: 3.0, .type: Option::Put, .strike: 90, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.30},
262 { .barrierType: Barrier::DownOut, .barrier: 95.0, .rebate: 3.0, .type: Option::Put, .strike: 100, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.30},
263 { .barrierType: Barrier::DownOut, .barrier: 95.0, .rebate: 3.0, .type: Option::Put, .strike: 110, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.30},
264 { .barrierType: Barrier::DownOut, .barrier: 100.0, .rebate: 3.0, .type: Option::Put, .strike: 90, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.30},
265 { .barrierType: Barrier::DownOut, .barrier: 100.0, .rebate: 3.0, .type: Option::Put, .strike: 100, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.30},
266 { .barrierType: Barrier::DownOut, .barrier: 100.0, .rebate: 3.0, .type: Option::Put, .strike: 110, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.30},
267 { .barrierType: Barrier::UpOut, .barrier: 105.0, .rebate: 3.0, .type: Option::Put, .strike: 90, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.30},
268 { .barrierType: Barrier::UpOut, .barrier: 105.0, .rebate: 3.0, .type: Option::Put, .strike: 100, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.30},
269 { .barrierType: Barrier::UpOut, .barrier: 105.0, .rebate: 3.0, .type: Option::Put, .strike: 110, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.30},
270
271 { .barrierType: Barrier::DownIn, .barrier: 95.0, .rebate: 3.0, .type: Option::Put, .strike: 90, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.30},
272 { .barrierType: Barrier::DownIn, .barrier: 95.0, .rebate: 3.0, .type: Option::Put, .strike: 100, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.30},
273 { .barrierType: Barrier::DownIn, .barrier: 95.0, .rebate: 3.0, .type: Option::Put, .strike: 110, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.30},
274 { .barrierType: Barrier::DownIn, .barrier: 100.0, .rebate: 3.0, .type: Option::Put, .strike: 90, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.30},
275 { .barrierType: Barrier::DownIn, .barrier: 100.0, .rebate: 3.0, .type: Option::Put, .strike: 100, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.30},
276 { .barrierType: Barrier::DownIn, .barrier: 100.0, .rebate: 3.0, .type: Option::Put, .strike: 110, .s: 100.0, .q: 0.04, .r: 0.08, .t: 1.00, .v: 0.15},
277 { .barrierType: Barrier::UpIn, .barrier: 105.0, .rebate: 3.0, .type: Option::Put, .strike: 90, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.30},
278 { .barrierType: Barrier::UpIn, .barrier: 105.0, .rebate: 3.0, .type: Option::Put, .strike: 100, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.30},
279 { .barrierType: Barrier::UpIn, .barrier: 105.0, .rebate: 3.0, .type: Option::Put, .strike: 110, .s: 100.0, .q: 0.04, .r: 0.08, .t: 0.50, .v: 0.30}
280 };
281
282 const DayCounter dc = Actual365Fixed();
283 const Date todaysDate(28, March, 2004);
284 const Date exerciseDate(28, March, 2005);
285 Settings::instance().evaluationDate() = todaysDate;
286
287 Handle<Quote> spot(
288 ext::shared_ptr<Quote>(new SimpleQuote(0.0)));
289 ext::shared_ptr<SimpleQuote> qRate(new SimpleQuote(0.0));
290 Handle<YieldTermStructure> qTS(flatRate(forward: qRate, dc));
291 ext::shared_ptr<SimpleQuote> rRate(new SimpleQuote(0.0));
292 Handle<YieldTermStructure> rTS(flatRate(forward: rRate, dc));
293 ext::shared_ptr<SimpleQuote> vol(new SimpleQuote(0.0));
294 Handle<BlackVolTermStructure> volTS(flatVol(volatility: vol, dc));
295
296 ext::shared_ptr<BlackScholesMertonProcess> bsProcess(
297 new BlackScholesMertonProcess(spot, qTS, rTS, volTS));
298
299 ext::shared_ptr<PricingEngine> analyticEngine(
300 new AnalyticBarrierEngine(bsProcess));
301
302 for (auto& value : values) {
303 Date exDate = todaysDate + timeToDays(t: value.t, daysPerYear: 365);
304 ext::shared_ptr<Exercise> exercise(new EuropeanExercise(exDate));
305
306 ext::dynamic_pointer_cast<SimpleQuote>(r: spot.currentLink())->setValue(value.s);
307 qRate->setValue(value.q);
308 rRate->setValue(value.r);
309 vol->setValue(value.v);
310
311 ext::shared_ptr<StrikedTypePayoff> payoff(new PlainVanillaPayoff(value.type, value.strike));
312
313 BarrierOption barrierOption(value.barrierType, value.barrier, value.rebate, payoff,
314 exercise);
315
316 const Real v0 = vol->value()*vol->value();
317 ext::shared_ptr<HestonProcess> hestonProcess(
318 new HestonProcess(rTS, qTS, spot, v0, 1.0, v0, 0.005, 0.0));
319
320 barrierOption.setPricingEngine(ext::shared_ptr<PricingEngine>(
321 new FdHestonBarrierEngine(ext::make_shared<HestonModel>(
322 args&: hestonProcess), 200, 101, 3)));
323
324 const Real calculatedHE = barrierOption.NPV();
325
326 barrierOption.setPricingEngine(analyticEngine);
327 const Real expected = barrierOption.NPV();
328
329 const Real tol = 0.0025;
330 if (std::fabs(x: calculatedHE - expected)/expected > tol) {
331 BOOST_ERROR("Failed to reproduce expected Heston npv"
332 << "\n calculated: " << calculatedHE
333 << "\n expected: " << expected
334 << "\n tolerance: " << tol);
335 }
336 }
337}
338
339void FdHestonTest::testFdmHestonBarrier() {
340
341 BOOST_TEST_MESSAGE("Testing FDM with barrier option for Heston model vs "
342 "Black-Scholes model...");
343
344 Handle<Quote> s0(ext::shared_ptr<Quote>(new SimpleQuote(100.0)));
345
346 Handle<YieldTermStructure> rTS(flatRate(forward: 0.05, dc: Actual365Fixed()));
347 Handle<YieldTermStructure> qTS(flatRate(forward: 0.0 , dc: Actual365Fixed()));
348
349 ext::shared_ptr<HestonProcess> hestonProcess(
350 new HestonProcess(rTS, qTS, s0, 0.04, 2.5, 0.04, 0.66, -0.8));
351
352 Settings::instance().evaluationDate() = Date(28, March, 2004);
353 Date exerciseDate(28, March, 2005);
354
355 ext::shared_ptr<Exercise> exercise(new EuropeanExercise(exerciseDate));
356
357 ext::shared_ptr<StrikedTypePayoff> payoff(new
358 PlainVanillaPayoff(Option::Call, 100));
359
360 BarrierOption barrierOption(Barrier::UpOut, 135, 0.0, payoff, exercise);
361
362 barrierOption.setPricingEngine(ext::shared_ptr<PricingEngine>(
363 new FdHestonBarrierEngine(ext::make_shared<HestonModel>(
364 args&: hestonProcess), 50, 400, 100)));
365
366 const Real tol = 0.01;
367 const Real npvExpected = 9.1530;
368 const Real deltaExpected = 0.5218;
369 const Real gammaExpected = -0.0354;
370
371 if (std::fabs(x: barrierOption.NPV() - npvExpected) > tol) {
372 BOOST_ERROR("Failed to reproduce expected npv"
373 << "\n calculated: " << barrierOption.NPV()
374 << "\n expected: " << npvExpected
375 << "\n tolerance: " << tol);
376 }
377 if (std::fabs(x: barrierOption.delta() - deltaExpected) > tol) {
378 BOOST_ERROR("Failed to reproduce expected delta"
379 << "\n calculated: " << barrierOption.delta()
380 << "\n expected: " << deltaExpected
381 << "\n tolerance: " << tol);
382 }
383 if (std::fabs(x: barrierOption.gamma() - gammaExpected) > tol) {
384 BOOST_ERROR("Failed to reproduce expected gamma"
385 << "\n calculated: " << barrierOption.gamma()
386 << "\n expected: " << gammaExpected
387 << "\n tolerance: " << tol);
388 }
389}
390
391void FdHestonTest::testFdmHestonAmerican() {
392
393 BOOST_TEST_MESSAGE("Testing FDM with American option in Heston model...");
394
395 Handle<Quote> s0(ext::shared_ptr<Quote>(new SimpleQuote(100.0)));
396
397 Handle<YieldTermStructure> rTS(flatRate(forward: 0.05, dc: Actual365Fixed()));
398 Handle<YieldTermStructure> qTS(flatRate(forward: 0.0 , dc: Actual365Fixed()));
399
400 ext::shared_ptr<HestonProcess> hestonProcess(
401 new HestonProcess(rTS, qTS, s0, 0.04, 2.5, 0.04, 0.66, -0.8));
402
403 Settings::instance().evaluationDate() = Date(28, March, 2004);
404 Date exerciseDate(28, March, 2005);
405
406 ext::shared_ptr<Exercise> exercise(new AmericanExercise(exerciseDate));
407
408 ext::shared_ptr<StrikedTypePayoff> payoff(new
409 PlainVanillaPayoff(Option::Put, 100));
410
411 VanillaOption option(payoff, exercise);
412 ext::shared_ptr<PricingEngine> engine(
413 new FdHestonVanillaEngine(ext::make_shared<HestonModel>(
414 args&: hestonProcess), 200, 100, 50));
415 option.setPricingEngine(engine);
416
417 const Real tol = 0.01;
418 const Real npvExpected = 5.66032;
419 const Real deltaExpected = -0.30065;
420 const Real gammaExpected = 0.02202;
421
422 if (std::fabs(x: option.NPV() - npvExpected) > tol) {
423 BOOST_ERROR("Failed to reproduce expected npv"
424 << "\n calculated: " << option.NPV()
425 << "\n expected: " << npvExpected
426 << "\n tolerance: " << tol);
427 }
428 if (std::fabs(x: option.delta() - deltaExpected) > tol) {
429 BOOST_ERROR("Failed to reproduce expected delta"
430 << "\n calculated: " << option.delta()
431 << "\n expected: " << deltaExpected
432 << "\n tolerance: " << tol);
433 }
434 if (std::fabs(x: option.gamma() - gammaExpected) > tol) {
435 BOOST_ERROR("Failed to reproduce expected gamma"
436 << "\n calculated: " << option.gamma()
437 << "\n expected: " << gammaExpected
438 << "\n tolerance: " << tol);
439 }
440}
441
442
443void FdHestonTest::testFdmHestonIkonenToivanen() {
444
445 BOOST_TEST_MESSAGE("Testing FDM Heston for Ikonen and Toivanen tests...");
446
447 /* check prices of american puts as given in:
448 From Efficient numerical methods for pricing American options under
449 stochastic volatility, Samuli Ikonen, Jari Toivanen,
450 http://users.jyu.fi/~tene/papers/reportB12-05.pdf
451 */
452 Handle<YieldTermStructure> rTS(flatRate(forward: 0.10, dc: Actual360()));
453 Handle<YieldTermStructure> qTS(flatRate(forward: 0.0 , dc: Actual360()));
454
455 Settings::instance().evaluationDate() = Date(28, March, 2004);
456 Date exerciseDate(26, June, 2004);
457
458 ext::shared_ptr<Exercise> exercise(new AmericanExercise(exerciseDate));
459
460 ext::shared_ptr<StrikedTypePayoff> payoff(new
461 PlainVanillaPayoff(Option::Put, 10));
462
463 VanillaOption option(payoff, exercise);
464
465 Real strikes[] = { 8, 9, 10, 11, 12 };
466 Real expected[] = { 2.00000, 1.10763, 0.520038, 0.213681, 0.082046 };
467 const Real tol = 0.001;
468
469 for (Size i=0; i < LENGTH(strikes); ++i) {
470 Handle<Quote> s0(ext::shared_ptr<Quote>(new SimpleQuote(strikes[i])));
471 ext::shared_ptr<HestonProcess> hestonProcess(
472 new HestonProcess(rTS, qTS, s0, 0.0625, 5, 0.16, 0.9, 0.1));
473
474 ext::shared_ptr<PricingEngine> engine(
475 new FdHestonVanillaEngine(ext::make_shared<HestonModel>(
476 args&: hestonProcess), 100, 400));
477 option.setPricingEngine(engine);
478
479 Real calculated = option.NPV();
480 if (std::fabs(x: calculated - expected[i]) > tol) {
481 BOOST_ERROR("Failed to reproduce expected npv"
482 << "\n strike: " << strikes[i]
483 << "\n calculated: " << calculated
484 << "\n expected: " << expected[i]
485 << "\n tolerance: " << tol);
486 }
487 }
488}
489
490void FdHestonTest::testFdmHestonBlackScholes() {
491
492 BOOST_TEST_MESSAGE("Testing FDM Heston with Black Scholes model...");
493
494 Settings::instance().evaluationDate() = Date(28, March, 2004);
495 Date exerciseDate(26, June, 2004);
496
497 Handle<YieldTermStructure> rTS(flatRate(forward: 0.10, dc: Actual360()));
498 Handle<YieldTermStructure> qTS(flatRate(forward: 0.0 , dc: Actual360()));
499 Handle<BlackVolTermStructure> volTS(
500 flatVol(today: rTS->referenceDate(), volatility: 0.25, dc: rTS->dayCounter()));
501
502 ext::shared_ptr<Exercise> exercise(new EuropeanExercise(exerciseDate));
503
504 ext::shared_ptr<StrikedTypePayoff> payoff(new
505 PlainVanillaPayoff(Option::Put, 10));
506
507 VanillaOption option(payoff, exercise);
508
509 Real strikes[] = { 8, 9, 10, 11, 12 };
510 const Real tol = 0.0001;
511
512 for (Real& strike : strikes) {
513 Handle<Quote> s0(ext::shared_ptr<Quote>(new SimpleQuote(strike)));
514
515 ext::shared_ptr<GeneralizedBlackScholesProcess> bsProcess(
516 new GeneralizedBlackScholesProcess(s0, qTS, rTS, volTS));
517
518 option.setPricingEngine(ext::shared_ptr<PricingEngine>(
519 new AnalyticEuropeanEngine(bsProcess)));
520
521 const Real expected = option.NPV();
522
523 ext::shared_ptr<HestonProcess> hestonProcess(
524 new HestonProcess(rTS, qTS, s0, 0.0625, 1, 0.0625, 0.0001, 0.0));
525
526 // Hundsdorfer scheme
527 option.setPricingEngine(ext::shared_ptr<PricingEngine>(
528 new FdHestonVanillaEngine(ext::make_shared<HestonModel>(
529 args&: hestonProcess),
530 100, 400, 3)));
531
532 Real calculated = option.NPV();
533 if (std::fabs(x: calculated - expected) > tol) {
534 BOOST_ERROR("Failed to reproduce expected npv"
535 << "\n strike: " << strike << "\n calculated: " << calculated
536 << "\n expected: " << expected << "\n tolerance: " << tol);
537 }
538
539 // Explicit scheme
540 option.setPricingEngine(ext::shared_ptr<PricingEngine>(
541 new FdHestonVanillaEngine(ext::make_shared<HestonModel>(
542 args&: hestonProcess),
543 4000, 400, 3, 0,
544 FdmSchemeDesc::ExplicitEuler())));
545
546 calculated = option.NPV();
547 if (std::fabs(x: calculated - expected) > tol) {
548 BOOST_ERROR("Failed to reproduce expected npv"
549 << "\n strike: " << strike << "\n calculated: " << calculated
550 << "\n expected: " << expected << "\n tolerance: " << tol);
551 }
552 }
553}
554
555
556
557void FdHestonTest::testFdmHestonEuropeanWithDividends() {
558
559 BOOST_TEST_MESSAGE("Testing FDM with European option with dividends in Heston model...");
560
561 Handle<Quote> s0(ext::shared_ptr<Quote>(new SimpleQuote(100.0)));
562
563 Handle<YieldTermStructure> rTS(flatRate(forward: 0.05, dc: Actual365Fixed()));
564 Handle<YieldTermStructure> qTS(flatRate(forward: 0.0 , dc: Actual365Fixed()));
565
566 auto hestonProcess = ext::make_shared<HestonProcess>(args&: rTS, args&: qTS, args&: s0, args: 0.04, args: 2.5, args: 0.04, args: 0.66, args: -0.8);
567
568 Settings::instance().evaluationDate() = Date(28, March, 2004);
569 Date exerciseDate(28, March, 2005);
570
571 auto exercise = ext::make_shared<AmericanExercise>(args&: exerciseDate);
572 auto payoff = ext::make_shared<PlainVanillaPayoff>(args: Option::Put, args: 100);
573
574 const std::vector<Real> dividends(1, 5);
575 const std::vector<Date> dividendDates(1, Date(28, September, 2004));
576
577 QL_DEPRECATED_DISABLE_WARNING
578 DividendVanillaOption option1(payoff, exercise, dividendDates, dividends);
579 QL_DEPRECATED_ENABLE_WARNING
580 ext::shared_ptr<PricingEngine> engine1(
581 new FdHestonVanillaEngine(ext::make_shared<HestonModel>(
582 args&: hestonProcess), 50, 100, 50));
583 option1.setPricingEngine(engine1);
584
585 const Real tol = 0.01;
586 const Real gammaTol = 0.001;
587 const Real npvExpected = 7.38216;
588 const Real deltaExpected = -0.397902;
589 const Real gammaExpected = 0.027747;
590
591 if (std::fabs(x: option1.NPV() - npvExpected) > tol) {
592 BOOST_ERROR("Failed to reproduce expected npv"
593 << "\n calculated: " << option1.NPV()
594 << "\n expected: " << npvExpected
595 << "\n tolerance: " << tol);
596 }
597 if (std::fabs(x: option1.delta() - deltaExpected) > tol) {
598 BOOST_ERROR("Failed to reproduce expected delta"
599 << "\n calculated: " << option1.delta()
600 << "\n expected: " << deltaExpected
601 << "\n tolerance: " << tol);
602 }
603 if (std::fabs(x: option1.gamma() - gammaExpected) > gammaTol) {
604 BOOST_ERROR("Failed to reproduce expected gamma"
605 << "\n calculated: " << option1.gamma()
606 << "\n expected: " << gammaExpected
607 << "\n tolerance: " << tol);
608 }
609
610
611 VanillaOption option2(payoff, exercise);
612 auto engine2 = ext::make_shared<FdHestonVanillaEngine>(
613 args: ext::make_shared<HestonModel>(args&: hestonProcess),
614 args: DividendVector(dividendDates, dividends),
615 args: 50, args: 100, args: 50);
616 option2.setPricingEngine(engine2);
617
618 if (std::fabs(x: option2.NPV() - npvExpected) > tol) {
619 BOOST_ERROR("Failed to reproduce expected npv"
620 << "\n calculated: " << option2.NPV()
621 << "\n expected: " << npvExpected
622 << "\n tolerance: " << tol);
623 }
624 if (std::fabs(x: option2.delta() - deltaExpected) > tol) {
625 BOOST_ERROR("Failed to reproduce expected delta"
626 << "\n calculated: " << option2.delta()
627 << "\n expected: " << deltaExpected
628 << "\n tolerance: " << tol);
629 }
630 if (std::fabs(x: option2.gamma() - gammaExpected) > gammaTol) {
631 BOOST_ERROR("Failed to reproduce expected gamma"
632 << "\n calculated: " << option2.gamma()
633 << "\n expected: " << gammaExpected
634 << "\n tolerance: " << tol);
635 }
636}
637
638namespace {
639 struct HestonTestData {
640 Real kappa;
641 Real theta;
642 Real sigma;
643 Real rho;
644 Real r;
645 Real q;
646 Real T;
647 Real K;
648 };
649}
650
651void FdHestonTest::testFdmHestonConvergence() {
652
653 /* convergence tests based on
654 ADI finite difference schemes for option pricing in the
655 Heston model with correlation, K.J. in t'Hout and S. Foulon
656 */
657
658 BOOST_TEST_MESSAGE("Testing FDM Heston convergence...");
659
660 HestonTestData values[] = {
661 { .kappa: 1.5 , .theta: 0.04 , .sigma: 0.3 , .rho: -0.9 , .r: 0.025 , .q: 0.0 , .T: 1.0 , .K: 100 },
662 { .kappa: 3.0 , .theta: 0.12 , .sigma: 0.04 , .rho: 0.6 , .r: 0.01 , .q: 0.04 , .T: 1.0 , .K: 100 },
663 { .kappa: 0.6067, .theta: 0.0707, .sigma: 0.2928, .rho: -0.7571, .r: 0.03 , .q: 0.0 , .T: 3.0 , .K: 100 },
664 { .kappa: 2.5 , .theta: 0.06 , .sigma: 0.5 , .rho: -0.1 , .r: 0.0507, .q: 0.0469, .T: 0.25, .K: 100 }
665 };
666
667 FdmSchemeDesc schemes[] = {
668 FdmSchemeDesc::Hundsdorfer(),
669 FdmSchemeDesc::ModifiedCraigSneyd(),
670 FdmSchemeDesc::ModifiedHundsdorfer(),
671 FdmSchemeDesc::CraigSneyd(),
672 FdmSchemeDesc::TrBDF2(),
673 FdmSchemeDesc::CrankNicolson(),
674 };
675
676 Size tn[] = { 60 };
677 Real v0[] = { 0.04 };
678
679 const Date todaysDate(28, March, 2004);
680 Settings::instance().evaluationDate() = todaysDate;
681
682 Handle<Quote> s0(ext::shared_ptr<Quote>(new SimpleQuote(75.0)));
683
684 for (const auto& scheme : schemes) {
685 for (auto& value : values) {
686 for (unsigned long j : tn) {
687 for (Real k : v0) {
688 Handle<YieldTermStructure> rTS(flatRate(forward: value.r, dc: Actual365Fixed()));
689 Handle<YieldTermStructure> qTS(flatRate(forward: value.q, dc: Actual365Fixed()));
690
691 ext::shared_ptr<HestonProcess> hestonProcess(new HestonProcess(
692 rTS, qTS, s0, k, value.kappa, value.theta, value.sigma, value.rho));
693
694 Date exerciseDate =
695 todaysDate + Period(static_cast<Integer>(value.T * 365), Days);
696 ext::shared_ptr<Exercise> exercise(
697 new EuropeanExercise(exerciseDate));
698
699 ext::shared_ptr<StrikedTypePayoff> payoff(
700 new PlainVanillaPayoff(Option::Call, value.K));
701
702 VanillaOption option(payoff, exercise);
703 ext::shared_ptr<PricingEngine> engine(new FdHestonVanillaEngine(
704 ext::make_shared<HestonModel>(args&: hestonProcess), j, 101, 51, 0, scheme));
705 option.setPricingEngine(engine);
706
707 const Real calculated = option.NPV();
708
709 ext::shared_ptr<PricingEngine> analyticEngine(
710 new AnalyticHestonEngine(
711 ext::make_shared<HestonModel>(
712 args&: hestonProcess), 144));
713
714 option.setPricingEngine(analyticEngine);
715 const Real expected = option.NPV();
716 if ( std::fabs(x: expected - calculated)/expected > 0.02
717 && std::fabs(x: expected - calculated) > 0.002) {
718 BOOST_ERROR("Failed to reproduce expected npv"
719 << "\n calculated: " << calculated
720 << "\n expected: " << expected
721 << "\n tolerance: " << 0.01);
722 }
723 }
724 }
725 }
726 }
727}
728
729void FdHestonTest::testFdmHestonIntradayPricing() {
730#ifdef QL_HIGH_RESOLUTION_DATE
731
732 BOOST_TEST_MESSAGE("Testing FDM Heston intraday pricing...");
733
734 const Option::Type type(Option::Put);
735 const Real underlying = 36;
736 const Real strike = underlying;
737 const Spread dividendYield = 0.00;
738 const Rate riskFreeRate = 0.06;
739 const Real v0 = 0.2;
740 const Real kappa = 1.0;
741 const Real theta = v0;
742 const Real sigma = 0.0065;
743 const Real rho = -0.75;
744 const DayCounter dayCounter = Actual365Fixed();
745
746 const Date maturity(17, May, 2014, 17, 30, 0);
747
748 const ext::shared_ptr<Exercise> europeanExercise(
749 new EuropeanExercise(maturity));
750 const ext::shared_ptr<StrikedTypePayoff> payoff(
751 new PlainVanillaPayoff(type, strike));
752 VanillaOption option(payoff, europeanExercise);
753
754 const Handle<Quote> s0(
755 ext::shared_ptr<Quote>(new SimpleQuote(underlying)));
756 RelinkableHandle<BlackVolTermStructure> flatVolTS;
757 RelinkableHandle<YieldTermStructure> flatTermStructure, flatDividendTS;
758 const ext::shared_ptr<HestonProcess> process(
759 new HestonProcess(flatTermStructure, flatDividendTS, s0,
760 v0, kappa, theta, sigma, rho));
761 const ext::shared_ptr<HestonModel> model(new HestonModel(process));
762 const ext::shared_ptr<PricingEngine> fdm(
763 new FdHestonVanillaEngine(model, 20, 100, 26, 0));
764 option.setPricingEngine(fdm);
765
766 const Real gammaExpected[] = {
767 1.46757, 1.54696, 1.6408, 1.75409, 1.89464,
768 2.07548, 2.32046, 2.67944, 3.28164, 4.64096 };
769
770 for (Size i = 0; i < 10; ++i) {
771 const Date now(17, May, 2014, 15, i*15, 0);
772 Settings::instance().evaluationDate() = now;
773
774 flatTermStructure.linkTo(ext::shared_ptr<YieldTermStructure>(
775 new FlatForward(now, riskFreeRate, dayCounter)));
776 flatDividendTS.linkTo(ext::shared_ptr<YieldTermStructure>(
777 new FlatForward(now, dividendYield, dayCounter)));
778
779 const Real gammaCalculated = option.gamma();
780 if (std::fabs(gammaCalculated - gammaExpected[i]) > 1e-4) {
781 BOOST_ERROR("unable to reproduce intraday gamma values at time "
782 << "\n timestamp : " << io::iso_datetime(now)
783 << "\n expiry : " << io::iso_datetime(maturity)
784 << "\n expected : " << gammaExpected[i]
785 << "\n calculated: "<< gammaCalculated);
786 }
787 }
788#endif
789}
790
791void FdHestonTest::testMethodOfLinesAndCN() {
792 BOOST_TEST_MESSAGE("Testing method of lines to solve Heston PDEs...");
793
794 const DayCounter dc = Actual365Fixed();
795 const Date today = Date(21, February, 2018);
796
797 Settings::instance().evaluationDate() = today;
798
799 const Handle<Quote> spot(ext::make_shared<SimpleQuote>(args: 100.0));
800 const Handle<YieldTermStructure> qTS(flatRate(today, forward: 0.0, dc));
801 const Handle<YieldTermStructure> rTS(flatRate(today, forward: 0.0, dc));
802
803 const Real v0 = 0.09;
804 const Real kappa = 1.0;
805 const Real theta = v0;
806 const Real sigma = 0.4;
807 const Real rho = -0.75;
808
809 const Date maturity = today + Period(3, Months);
810
811 const ext::shared_ptr<HestonModel> model(
812 ext::make_shared<HestonModel>(
813 args: ext::make_shared<HestonProcess>(
814 args: rTS, args: qTS, args: spot, args: v0, args: kappa, args: theta, args: sigma, args: rho)));
815
816 const Size xGrid = 21;
817 const Size vGrid = 7;
818
819 const ext::shared_ptr<PricingEngine> fdmDefault(
820 ext::make_shared<FdHestonVanillaEngine>(args: model, args: 10, args: xGrid, args: vGrid, args: 0));
821
822 const ext::shared_ptr<PricingEngine> fdmMol(
823 ext::make_shared<FdHestonVanillaEngine>(
824 args: model, args: 10, args: xGrid, args: vGrid, args: 0, args: FdmSchemeDesc::MethodOfLines()));
825
826 const ext::shared_ptr<PlainVanillaPayoff> payoff =
827 ext::make_shared<PlainVanillaPayoff>(args: Option::Put, args: spot->value());
828
829 VanillaOption option(
830 payoff, ext::make_shared<AmericanExercise>(args: maturity));
831
832 option.setPricingEngine(fdmMol);
833 const Real calculatedMoL = option.NPV();
834
835 option.setPricingEngine(fdmDefault);
836 const Real expected = option.NPV();
837
838 const Real tol = 0.005;
839 const Real diffMoL = std::fabs(x: expected - calculatedMoL);
840
841 if (diffMoL > tol) {
842 BOOST_FAIL("Failed to reproduce european option values with MOL"
843 << "\n calculated: " << calculatedMoL
844 << "\n expected: " << expected
845 << "\n difference: " << diffMoL
846 << "\n tolerance: " << tol);
847 }
848
849 const ext::shared_ptr<PricingEngine> fdmCN(
850 ext::make_shared<FdHestonVanillaEngine>(
851 args: model, args: 10, args: xGrid, args: vGrid, args: 0, args: FdmSchemeDesc::CrankNicolson()));
852 option.setPricingEngine(fdmCN);
853
854 const Real calculatedCN = option.NPV();
855 const Real diffCN = std::fabs(x: expected - calculatedCN);
856
857 if (diffCN > tol) {
858 BOOST_FAIL("Failed to reproduce european option values with Crank-Nicolson"
859 << "\n calculated: " << calculatedCN
860 << "\n expected: " << expected
861 << "\n difference: " << diffCN
862 << "\n tolerance: " << tol);
863 }
864
865 BarrierOption barrierOption(
866 Barrier::DownOut, 85.0, 10.0,
867 payoff, ext::make_shared<EuropeanExercise>(args: maturity));
868
869 barrierOption.setPricingEngine(
870 ext::make_shared<FdHestonBarrierEngine>(args: model, args: 100, args: 31, args: 11));
871
872 const Real expectedBarrier = barrierOption.NPV();
873
874 barrierOption.setPricingEngine(
875 ext::make_shared<FdHestonBarrierEngine>(args: model, args: 100, args: 31, args: 11, args: 0,
876 args: FdmSchemeDesc::MethodOfLines()));
877
878 const Real calculatedBarrierMoL = barrierOption.NPV();
879
880 const Real barrierTol = 0.01;
881 const Real barrierDiffMoL = std::fabs(x: expectedBarrier - calculatedBarrierMoL);
882
883 if (barrierDiffMoL > barrierTol) {
884 BOOST_FAIL("Failed to reproduce barrier option values with MOL"
885 << "\n calculated: " << calculatedBarrierMoL
886 << "\n expected: " << expectedBarrier
887 << "\n difference: " << barrierDiffMoL
888 << "\n tolerance: " << barrierTol);
889 }
890
891 barrierOption.setPricingEngine(
892 ext::make_shared<FdHestonBarrierEngine>(args: model, args: 100, args: 31, args: 11, args: 0,
893 args: FdmSchemeDesc::CrankNicolson()));
894
895 const Real calculatedBarrierCN = barrierOption.NPV();
896 const Real barrierDiffCN = std::fabs(x: expectedBarrier - calculatedBarrierCN);
897
898 if (barrierDiffCN > barrierTol) {
899 BOOST_FAIL("Failed to reproduce barrier option values with Crank-Nicolson"
900 << "\n calculated: " << calculatedBarrierCN
901 << "\n expected: " << expectedBarrier
902 << "\n difference: " << barrierDiffCN
903 << "\n tolerance: " << barrierTol);
904 }
905}
906
907void FdHestonTest::testSpuriousOscillations() {
908 BOOST_TEST_MESSAGE("Testing for spurious oscillations when "
909 "solving the Heston PDEs...");
910
911 const DayCounter dc = Actual365Fixed();
912 const Date today = Date(7, June, 2018);
913
914 Settings::instance().evaluationDate() = today;
915
916 const Handle<Quote> spot(ext::make_shared<SimpleQuote>(args: 100.0));
917 const Handle<YieldTermStructure> qTS(flatRate(today, forward: 0.1, dc));
918 const Handle<YieldTermStructure> rTS(flatRate(today, forward: 0.0, dc));
919
920 const Real v0 = 0.005;
921 const Real kappa = 1.0;
922 const Real theta = 0.005;
923 const Real sigma = 0.4;
924 const Real rho = -0.75;
925
926 const Date maturity = today + Period(1, Years);
927
928 const ext::shared_ptr<HestonProcess> process =
929 ext::make_shared<HestonProcess>(
930 args: rTS, args: qTS, args: spot, args: v0, args: kappa, args: theta, args: sigma, args: rho);
931
932 const ext::shared_ptr<HestonModel> model =
933 ext::make_shared<HestonModel>(args: process);
934
935 const ext::shared_ptr<FdHestonVanillaEngine> hestonEngine(
936 ext::make_shared<FdHestonVanillaEngine>(
937 args: model, args: 6, args: 200, args: 13, args: 0, args: FdmSchemeDesc::TrBDF2()));
938
939 VanillaOption option(
940 ext::make_shared<PlainVanillaPayoff>(args: Option::Call, args: spot->value()),
941 ext::make_shared<EuropeanExercise>(args: maturity));
942
943 option.setupArguments(hestonEngine->getArguments());
944
945 const ext::tuple<FdmSchemeDesc, std::string, bool> descs[] = {
946 ext::make_tuple(args: FdmSchemeDesc::CraigSneyd(), args: "Craig-Sneyd", args: true),
947 ext::make_tuple(args: FdmSchemeDesc::Hundsdorfer(), args: "Hundsdorfer", args: true),
948 ext::make_tuple(
949 args: FdmSchemeDesc::ModifiedHundsdorfer(), args: "Mod. Hundsdorfer", args: true),
950 ext::make_tuple(args: FdmSchemeDesc::Douglas(), args: "Douglas", args: true),
951 ext::make_tuple(args: FdmSchemeDesc::CrankNicolson(), args: "Crank-Nicolson", args: true),
952 ext::make_tuple(args: FdmSchemeDesc::ImplicitEuler(), args: "Implicit", args: false),
953 ext::make_tuple(args: FdmSchemeDesc::TrBDF2(), args: "TR-BDF2", args: false)
954 };
955
956 for (const auto& desc : descs) {
957 const ext::shared_ptr<FdmHestonSolver> solver = ext::make_shared<FdmHestonSolver>(
958 args: Handle<HestonProcess>(process), args: hestonEngine->getSolverDesc(equityScaleFactor: 1.0), args: ext::get<0>(t: desc));
959
960 std::vector<Real> gammas;
961 for (Real x=99; x < 101.001; x+=0.1) {
962 gammas.push_back(x: solver->gammaAt(s: x, v: v0));
963 }
964
965 Real maximum = QL_MIN_REAL;
966 for (Size i=1; i < gammas.size(); ++i) {
967 const Real diff = std::fabs(x: gammas[i] - gammas[i-1]);
968 if (diff > maximum)
969 maximum = diff;
970 }
971
972 const Real tol = 0.01;
973 const bool hasSpuriousOscillations = maximum > tol;
974
975 if (hasSpuriousOscillations != ext::get<2>(t: desc)) {
976 BOOST_ERROR("unable to reproduce spurious oscillation behaviour "
977 << "\n scheme name : " << ext::get<1>(desc)
978 << "\n oscillations observed: " << hasSpuriousOscillations
979 << "\n oscillations expected: " << ext::get<2>(desc));
980 }
981 }
982}
983
984
985void FdHestonTest::testAmericanCallPutParity() {
986 BOOST_TEST_MESSAGE("Testing call/put parity for American options "
987 "under the Heston model...");
988
989 // A. Battauz, M. De Donno,m A. Sbuelz:
990 // The put-call symmetry for American options in
991 // the Heston stochastic volatility model
992
993 const DayCounter dc = Actual365Fixed();
994 const Date today = Date(15, April, 2022);
995
996 Settings::instance().evaluationDate() = today;
997
998 struct OptionSpec {
999 Real spot;
1000 Real strike;
1001 Size maturityInDays;
1002 Real r, q;
1003 Real v0, kappa, theta, sig, rho;
1004 };
1005
1006 auto buildStochProcess = [&dc](const OptionSpec& testCase) {
1007 return ext::make_shared<HestonProcess>(
1008 args: Handle<YieldTermStructure>(flatRate(forward: testCase.r, dc)),
1009 args: Handle<YieldTermStructure>(flatRate(forward: testCase.q, dc)),
1010 args: Handle<Quote>(ext::make_shared<SimpleQuote>(args: testCase.spot)),
1011 args: testCase.v0, args: testCase.kappa,
1012 args: testCase.theta, args: testCase.sig, args: testCase.rho
1013 );
1014 };
1015
1016 const OptionSpec testCaseSpecs[] = {
1017 {.spot: 100.0, .strike: 90.0, .maturityInDays: 365, .r: 0.02, .q: 0.15, .v0: 0.25, .kappa: 1.0, .theta: 0.09, .sig: 0.5, .rho: -0.75},
1018 {.spot: 100.0, .strike: 90.0, .maturityInDays: 365, .r: 0.05, .q: 0.20, .v0: 0.5, .kappa: 1.0, .theta: 0.05, .sig: 0.75, .rho: -0.9}
1019 };
1020
1021 const Size xGrid = 200;
1022 const Size vGrid = 25;
1023 const Size timeStepsPerYear = 50;
1024
1025 for (const auto& testCaseSpec: testCaseSpecs) {
1026 const auto maturityDate =
1027 today + Period(testCaseSpec.maturityInDays, Days);
1028 const Time maturityTime = dc.yearFraction(d1: today, d2: maturityDate);
1029 const Size tGrid = Size(maturityTime * timeStepsPerYear);
1030
1031 const auto exercise =
1032 ext::make_shared<AmericanExercise>(args: today, args: maturityDate);
1033
1034 VanillaOption callOption(
1035 ext::make_shared<PlainVanillaPayoff>(
1036 args: Option::Call, args: testCaseSpec.strike),
1037 exercise
1038 );
1039
1040 callOption.setPricingEngine(
1041 ext::make_shared<FdHestonVanillaEngine>(
1042 args: ext::make_shared<HestonModel>(
1043 args: buildStochProcess(testCaseSpec)),
1044 args: tGrid, args: xGrid, args: vGrid
1045 )
1046 );
1047
1048 const Real callNpv = callOption.NPV();
1049
1050 OptionSpec putOptionSpec = {
1051 .spot: testCaseSpec.strike,
1052 .strike: testCaseSpec.spot,
1053 .maturityInDays: testCaseSpec.maturityInDays,
1054 .r: testCaseSpec.q,
1055 .q: testCaseSpec.r,
1056 .v0: testCaseSpec.v0,
1057 .kappa: testCaseSpec.kappa - testCaseSpec.sig*testCaseSpec.rho,
1058 .theta: testCaseSpec.kappa*testCaseSpec.theta/
1059 (testCaseSpec.kappa - testCaseSpec.sig*testCaseSpec.rho),
1060 .sig: testCaseSpec.sig,
1061 .rho: -testCaseSpec.rho
1062 };
1063
1064 VanillaOption putOption(
1065 ext::make_shared<PlainVanillaPayoff>(
1066 args: Option::Put, args&: putOptionSpec.strike),
1067 exercise
1068 );
1069
1070 putOption.setPricingEngine(
1071 ext::make_shared<FdHestonVanillaEngine>(
1072 args: ext::make_shared<HestonModel>(
1073 args: buildStochProcess(putOptionSpec)),
1074 args: tGrid, args: xGrid, args: vGrid
1075 )
1076 );
1077
1078 const Real putNpv = putOption.NPV();
1079
1080 const Real diff = std::fabs(x: putNpv -callNpv);
1081 const Real tol = 0.025;
1082
1083 if (diff > tol) {
1084 BOOST_FAIL("failed to reproduce American call/put parity"
1085 << "\n Put NPV : " << putNpv
1086 << "\n Call NPV : " << callNpv
1087 << "\n difference: " << diff
1088 << "\n tolerance : " << tol);
1089 }
1090 }
1091}
1092
1093test_suite* FdHestonTest::suite(SpeedLevel speed) {
1094 auto* suite = BOOST_TEST_SUITE("Finite Difference Heston tests");
1095
1096 suite->add(QUANTLIB_TEST_CASE(&FdHestonTest::testFdmHestonVarianceMesher));
1097 suite->add(QUANTLIB_TEST_CASE(&FdHestonTest::testFdmHestonBarrier));
1098 suite->add(QUANTLIB_TEST_CASE(&FdHestonTest::testFdmHestonAmerican));
1099 suite->add(QUANTLIB_TEST_CASE(&FdHestonTest::testFdmHestonIkonenToivanen));
1100 suite->add(QUANTLIB_TEST_CASE(&FdHestonTest::testFdmHestonEuropeanWithDividends));
1101 #ifdef QL_HIGH_RESOLUTION_DATE
1102 suite->add(QUANTLIB_TEST_CASE(&FdHestonTest::testFdmHestonIntradayPricing));
1103 #endif
1104 suite->add(QUANTLIB_TEST_CASE(&FdHestonTest::testMethodOfLinesAndCN));
1105 suite->add(QUANTLIB_TEST_CASE(&FdHestonTest::testSpuriousOscillations));
1106 suite->add(QUANTLIB_TEST_CASE(&FdHestonTest::testAmericanCallPutParity));
1107 suite->add(QUANTLIB_TEST_CASE(&FdHestonTest::testFdmHestonBlackScholes));
1108
1109 if (speed <= Fast) {
1110 suite->add(QUANTLIB_TEST_CASE(&FdHestonTest::testFdmHestonConvergence));
1111 suite->add(QUANTLIB_TEST_CASE(&FdHestonTest::testFdmHestonBarrierVsBlackScholes));
1112 }
1113
1114 return suite;
1115}
1116
1117

source code of quantlib/test-suite/fdheston.cpp