1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2020 Jack Gillett
5
6 This file is part of QuantLib, a free-software/open-source library
7 for financial quantitative analysts and developers - http://quantlib.org/
8
9 QuantLib is free software: you can redistribute it and/or modify it
10 under the terms of the QuantLib license. You should have received a
11 copy of the license along with this program; if not, please email
12 <quantlib-dev@lists.sf.net>. The license is also available online at
13 <http://quantlib.org/license.shtml>.
14
15 This program is distributed in the hope that it will be useful, but WITHOUT
16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17 FOR A PARTICULAR PURPOSE. See the license for more details.
18*/
19
20#include <ql/experimental/forward/analytichestonforwardeuropeanengine.hpp>
21#include <complex>
22#include <utility>
23
24namespace QuantLib {
25
26
27 class P12Integrand {
28 private:
29 ext::shared_ptr<AnalyticHestonEngine>& engine_;
30 Real logK_, phiRightLimit_;
31 Time tenor_;
32 std::complex<Real> i_, adj_;
33 public:
34 P12Integrand(ext::shared_ptr<AnalyticHestonEngine>& engine,
35 Real logK,
36 Time tenor,
37 bool P1, // true for P1, false for P2
38 Real phiRightLimit = 100) : engine_(engine), logK_(logK),
39 phiRightLimit_(phiRightLimit), tenor_(tenor), i_(std::complex<Real>(0.0, 1.0)) {
40
41 // Only difference between P1 and P2 integral is the additional term in the chF evaluation
42 if (P1) {
43 adj_ = std::complex<Real>(0.0, -1.0);
44 } else {
45 adj_ = std::complex<Real>(0.0, 0.0);
46 }
47 }
48
49 // QL Gaussian Quadrature - map phi from [-1 to 1] to {0, phiRightLimit]
50 Real operator()(Real phi) const {
51 Real phiDash = (0.5+1e-8+0.5*phi) * phiRightLimit_; // Map phi to full range
52 return 0.5*phiRightLimit_*std::real(z: (std::exp(z: -phiDash*logK_*i_) / (phiDash*i_)) * engine_->chF(z: phiDash+adj_, t: tenor_));
53 }
54 };
55
56
57 class P12HatIntegrand {
58 private:
59 Time tenor_, resetTime_;
60 Handle<Quote>& s0_;
61 bool P1_;
62 Real logK_, phiRightLimit_, nuRightLimit_;
63 const AnalyticHestonForwardEuropeanEngine* const parent_;
64 GaussLegendreIntegration innerIntegrator_;
65 public:
66 P12HatIntegrand(Time tenor,
67 Time resetTime,
68 Handle<Quote>& s0,
69 Real logK,
70 bool P1, // true for P1, false for P2
71 const AnalyticHestonForwardEuropeanEngine* const parent,
72 Real phiRightLimit,
73 Real nuRightLimit) : tenor_(tenor), resetTime_(resetTime),
74 s0_(s0), P1_(P1), logK_(logK), phiRightLimit_(phiRightLimit),
75 nuRightLimit_(nuRightLimit), parent_(parent), innerIntegrator_(128) {}
76 Real operator()(Real nu) const {
77
78 // Rescale nu to [-1, 1]
79 Real nuDash = nuRightLimit_ * (0.5 * nu + 0.5 + 1e-8);
80
81 // Calculate the chF from var(t) to expiry
82 ext::shared_ptr<AnalyticHestonEngine> engine = parent_->forwardChF(spotReset&: s0_, varReset: nuDash);
83 P12Integrand pIntegrand(engine, logK_, tenor_, P1_, phiRightLimit_);
84 Real p1Integral = innerIntegrator_(pIntegrand);
85
86 // Calculate the value of the propagator to nu
87 Real propagator = parent_->propagator(resetTime: resetTime_, varReset: nuDash);
88
89 // Take the product, and scale integral's value back up to [0, right_lim]
90 return propagator * (0.5 + p1Integral/M_PI);
91 }
92 };
93
94
95 AnalyticHestonForwardEuropeanEngine::AnalyticHestonForwardEuropeanEngine(
96 ext::shared_ptr<HestonProcess> process, Size integrationOrder)
97 : process_(std::move(process)), integrationOrder_(integrationOrder), outerIntegrator_(128) {
98
99 v0_ = process_->v0();
100 rho_ = process_->rho();
101 kappa_ = process_->kappa();
102 theta_ = process_->theta();
103 sigma_ = process_->sigma();
104 s0_ = process_->s0();
105
106 QL_REQUIRE(sigma_ > 0.1,
107 "Very low values (<~10%) for Heston Vol-of-Vol cause numerical issues" \
108 "in this implementation of the propagator function, try using" \
109 "MCForwardEuropeanHestonEngine Monte-Carlo engine instead");
110
111 riskFreeRate_ = process_->riskFreeRate();
112 dividendYield_ = process_->dividendYield();
113
114 // Some of the required constant intermediate variables can be calculated now
115 kappaHat_ = kappa_ - rho_ * sigma_;
116 thetaHat_ = kappa_ * theta_ / kappaHat_;
117 R_ = 4 * kappaHat_ * thetaHat_ / (sigma_ * sigma_);
118 }
119
120
121 void AnalyticHestonForwardEuropeanEngine::calculate() const {
122 // This is a european option pricer
123 QL_REQUIRE(this->arguments_.exercise->type() == Exercise::European,
124 "not an European option");
125
126 // We only price plain vanillas
127 ext::shared_ptr<PlainVanillaPayoff> payoff =
128 ext::dynamic_pointer_cast<PlainVanillaPayoff>(r: this->arguments_.payoff);
129 QL_REQUIRE(payoff, "non plain vanilla payoff given");
130
131 Time resetTime = this->process_->time(this->arguments_.resetDate);
132 Time expiryTime = this->process_->time(this->arguments_.exercise->lastDate());
133 Time tenor = expiryTime - resetTime;
134 Real moneyness = this->arguments_.moneyness;
135
136 // K needs to be scaled to forward AT RESET TIME, not spot...
137 Real expiryDcf = riskFreeRate_->discount(t: expiryTime);
138 Real resetDcf = riskFreeRate_->discount(t: resetTime);
139 Real expiryDividendDiscount = dividendYield_->discount(t: expiryTime);
140 Real resetDividendDiscount = dividendYield_->discount(t: resetTime);
141 Real expiryRatio = expiryDcf / expiryDividendDiscount;
142 Real resetRatio = resetDcf / resetDividendDiscount;
143
144 QL_REQUIRE(resetTime >= 0.0, "Reset Date cannot be in the past");
145 QL_REQUIRE(expiryTime >= 0.0, "Expiry Date cannot be in the past");
146
147 // Use some heuristics to decide upon phiRightLimit and nuRightLimit
148 Real phiRightLimit = 100.0;
149 Real nuRightLimit = std::max(a: 2.0, b: 10.0 * (1+std::max(a: 0.0, b: rho_)) * sigma_ * std::sqrt(x: resetTime * std::max(a: v0_, b: theta_)));
150
151 // do the 2D integral calculation. For very short times, we just fall back on the standard
152 // calculation, both for accuracy and because tStar==0 causes some numerical issues...
153 std::pair<Real, Real> P1HatP2Hat;
154 if (resetTime <= 1e-3) {
155 Handle<Quote> tempQuote(ext::shared_ptr<Quote>(new SimpleQuote(s0_->value())));
156 P1HatP2Hat = calculateP1P2(t: tenor, St&: tempQuote, K: moneyness * s0_->value(), ratio: expiryRatio, phiRightLimit);
157 } else {
158 P1HatP2Hat = calculateP1P2Hat(tenor, resetTime, K: moneyness, ratio: expiryRatio/resetRatio, phiRightLimit, nuRightLimit);
159 }
160
161 // Apply the payoff functions
162 Real value = 0.0;
163 Real F = s0_->value() / expiryRatio;
164 switch (payoff->optionType()){
165 case Option::Call:
166 value = expiryDcf * (F*P1HatP2Hat.first - moneyness*s0_->value()*P1HatP2Hat.second/resetRatio);
167 break;
168 case Option::Put:
169 value = expiryDcf * (moneyness*s0_->value()*(1-P1HatP2Hat.second)/resetRatio - F*(1-P1HatP2Hat.first));
170 break;
171 default:
172 QL_FAIL("unknown option type");
173 }
174
175 results_.value = value;
176
177 results_.additionalResults["dcf"] = expiryDcf;
178 results_.additionalResults["qf"] = expiryDividendDiscount;
179 results_.additionalResults["expiryRatio"] = expiryRatio;
180 results_.additionalResults["resetRatio"] = resetRatio;
181 results_.additionalResults["moneyness"] = moneyness;
182 results_.additionalResults["s0"] = s0_->value();
183 results_.additionalResults["fwd"] = F;
184 results_.additionalResults["resetTime"] = resetTime;
185 results_.additionalResults["expiryTime"] = expiryTime;
186 results_.additionalResults["P1Hat"] = P1HatP2Hat.first;
187 results_.additionalResults["P2Hat"] = P1HatP2Hat.second;
188 results_.additionalResults["phiRightLimit"] = phiRightLimit;
189 results_.additionalResults["nuRightLimit"] = nuRightLimit;
190 }
191
192
193 std::pair<Real, Real> AnalyticHestonForwardEuropeanEngine::calculateP1P2Hat(Time tenor,
194 Time resetTime,
195 Real moneyness,
196 Real ratio,
197 Real phiRightLimit,
198 Real nuRightLimit) const {
199
200 Handle<Quote> unitQuote(ext::shared_ptr<Quote>(new SimpleQuote(1.0)));
201
202 // Re-expressing moneyness in terms of the forward here (strike fixes to spot, but in
203 // our pricing calculation we need to compare it to the future at expiry)
204 Real logMoneyness = std::log(x: moneyness*ratio);
205
206 P12HatIntegrand p1HatIntegrand(tenor, resetTime, unitQuote, logMoneyness, true, this, phiRightLimit, nuRightLimit);
207 P12HatIntegrand p2HatIntegrand(tenor, resetTime, unitQuote, logMoneyness, false, this, phiRightLimit, nuRightLimit);
208
209 Real p1HatIntegral = 0.5 * nuRightLimit * outerIntegrator_(p1HatIntegrand);
210 Real p2HatIntegral = 0.5 * nuRightLimit * outerIntegrator_(p2HatIntegrand);
211
212 std::pair<Real, Real> P1HatP2Hat(p1HatIntegral, p2HatIntegral);
213
214 return P1HatP2Hat;
215 }
216
217
218 Real AnalyticHestonForwardEuropeanEngine::propagator(Time resetTime,
219 Real varReset) const {
220 Real B, Lambda, term1, term2, term3;
221
222 B = 4 * kappaHat_ / (sigma_ * sigma_ * (1 - std::exp(x: -kappaHat_ * resetTime)));
223 Lambda = B * std::exp(x: -kappaHat_ * resetTime) * v0_;
224
225 // Now construct equation (18) from the paper term-by-term
226 term1 = std::exp(x: -0.5*(B * varReset + Lambda)) * B / 2;
227 term2 = std::pow(x: B * varReset / Lambda, y: 0.5*(R_/2 - 1));
228 term3 = modifiedBesselFunction_i(nu: Real(R_/2 - 1),x: Real(std::sqrt(x: Lambda * B * varReset)));
229
230 return term1 * term2 * term3;
231 }
232
233 ext::shared_ptr<AnalyticHestonEngine> AnalyticHestonForwardEuropeanEngine::forwardChF(
234 Handle<Quote>& spotReset,
235 Real varReset) const {
236
237 // Probably a wasteful implementation here, could be improved by importing
238 // only the CF-generating parts of the AnalyticHestonEngine (currently private)
239 ext::shared_ptr<HestonProcess> hestonProcess(new
240 HestonProcess(riskFreeRate_, dividendYield_, spotReset,
241 varReset, kappa_, theta_, sigma_, rho_));
242
243 ext::shared_ptr<HestonModel> hestonModel(new HestonModel(hestonProcess));
244
245 ext::shared_ptr<AnalyticHestonEngine> analyticHestonEngine(
246 new AnalyticHestonEngine(hestonModel, integrationOrder_));
247
248 // Not sure how to pass only the chF, so just pass the whole thing for now!
249 return analyticHestonEngine;
250 }
251
252
253 std::pair<Real, Real> AnalyticHestonForwardEuropeanEngine::calculateP1P2(Time tenor,
254 Handle<Quote>& St,
255 Real K,
256 Real ratio,
257 Real phiRightLimit) const {
258
259 ext::shared_ptr<AnalyticHestonEngine> engine = forwardChF(spotReset&: St, varReset: v0_);
260 Real logK = std::log(x: K*ratio/St->value());
261
262 // Integrate the CF and the complex integrand over positive phi
263 GaussLegendreIntegration integrator = GaussLegendreIntegration(128);
264 P12Integrand p1Integrand(engine, logK, tenor, true, phiRightLimit);
265 P12Integrand p2Integrand(engine, logK, tenor, false, phiRightLimit);
266
267 Real p1Integral = integrator(p1Integrand);
268 Real p2Integral = integrator(p2Integrand);
269
270 std::pair<Real, Real> P1P2(0.5 + p1Integral/M_PI, 0.5 + p2Integral/M_PI);
271
272 return P1P2;
273 }
274}
275

source code of quantlib/ql/experimental/forward/analytichestonforwardeuropeanengine.cpp