QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
analytic_cont_geom_av_price_heston.cpp
Go to the documentation of this file.
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
21#include <utility>
22
23namespace QuantLib {
24
25 class AnalyticContinuousGeometricAveragePriceAsianHestonEngine::Integrand {
26 private:
27 Real t_ = 0.0, T_, K_, logK_;
28 Size cutoff_;
29 const AnalyticContinuousGeometricAveragePriceAsianHestonEngine* const parent_;
30 Real xiRightLimit_;
31 std::complex<Real> i_;
32
33 public:
34 Integrand(Real T,
35 Size cutoff,
36 Real K,
37 const AnalyticContinuousGeometricAveragePriceAsianHestonEngine* const parent,
38 Real xiRightLimit)
39 : T_(T), K_(K), logK_(std::log(K)), cutoff_(cutoff), parent_(parent),
40 xiRightLimit_(xiRightLimit), i_(std::complex<Real>(0.0, 1.0)) {}
41
42 Real operator()(Real xi) const {
43 Real xiDash = (0.5+1e-8+0.5*xi) * xiRightLimit_; // Map xi to full range
44
45 std::complex<Real> inner1 = parent_->Phi(1.0 + xiDash*i_, 0, T_, t_, cutoff_);
46 std::complex<Real> inner2 = - K_*parent_->Phi(xiDash*i_, 0, T_, t_, cutoff_);
47
48 return 0.5*xiRightLimit_*std::real((inner1 + inner2) * std::exp(-xiDash*logK_*i_) / (xiDash*i_));
49 }
50 };
51
52 class AnalyticContinuousGeometricAveragePriceAsianHestonEngine::DcfIntegrand {
53 private:
54 Real t_, T_, denominator_;
55 const Handle<YieldTermStructure> riskFreeRate_;
56 const Handle<YieldTermStructure> dividendYield_;
57 public:
58 DcfIntegrand(Real t,
59 Real T,
60 Handle<YieldTermStructure> riskFreeRate,
61 Handle<YieldTermStructure> dividendYield)
62 : t_(t), T_(T), riskFreeRate_(std::move(riskFreeRate)),
63 dividendYield_(std::move(dividendYield)) {
64 denominator_ = std::log(riskFreeRate_->discount(t_)) - std::log(dividendYield_->discount(t_));
65 }
66
67 Real operator()(Real u) const {
68 Real uDash = (0.5+1e-8+0.5*u) * (T_ - t_) + t_; // Map u to full range
69 return 0.5*(T_ - t_)*(-std::log(riskFreeRate_->discount(uDash))
70 + std::log(dividendYield_->discount(uDash)) + denominator_);
71 }
72 };
73
74
77 ext::shared_ptr<HestonProcess> process, Size summationCutoff, Real xiRightLimit)
78 : process_(std::move(process)), summationCutoff_(summationCutoff), xiRightLimit_(xiRightLimit),
79 integrator_(128) {
81
82 v0_ = process_->v0();
83 rho_ = process_->rho();
84 kappa_ = process_->kappa();
85 theta_ = process_->theta();
86 sigma_ = process_->sigma();
87 s0_ = process_->s0();
88
89 riskFreeRate_ = process_->riskFreeRate();
90 dividendYield_ = process_->dividendYield();
91
92 // Some of the required constant intermediate variables can be calculated now
93 // (although anything depending on T will need to be calculated dynamically later)
94 a1_ = 2.0 * v0_ / (sigma_*sigma_);
95 a2_ = 2.0 * kappa_ * theta_ / (sigma_*sigma_);
96 }
97
99 const std::complex<Real>& s, const std::complex<Real>& w, Real T) const {
100 return s*s*(1-rho_*rho_)/(2*T*T);
101 }
102
104 const std::complex<Real>& s, const std::complex<Real>& w, Real T) const {
105 return s*(2*rho_*kappa_ - sigma_)/(2*sigma_*T) + s*w*(1-rho_*rho_)/T;
106 }
107
109 const std::complex<Real>& s, const std::complex<Real>& w, Real T) const {
110 return s*rho_/(sigma_*T) + 0.5*w*(2*rho_*kappa_ - sigma_)/sigma_ + 0.5*w*w*(1-rho_*rho_);
111 }
112
114 const std::complex<Real>& s, const std::complex<Real>& w) const {
115 return w*rho_/sigma_;
116 }
117
118 std::complex<Real> AnalyticContinuousGeometricAveragePriceAsianHestonEngine::f(const std::complex<Real>& z1,
119 const std::complex<Real>& z2,
120 const std::complex<Real>& z3,
121 const std::complex<Real>& z4,
122 int n, // Can't use Size here as n can be negative
123 Real tau) const {;
124 std::complex<Real> result;
125
126 // This equation is highly recursive, use dynamic programming with a mutable variable
127 // to record the results of previous calls
128 if (n<2) {
129 if (n<0) {
130 result = 0.0;
131 } else if (n==0) {
132 result = 1.0;
133 } else {
134 result = 0.5*(kappa_ - z4*sigma_*sigma_)*tau;
135 }
136 } else {
137 std::complex<Real> fMinusN[4];
138 Real prefactor = -0.5*sigma_*sigma_*tau*tau / (n*(n-1));
139
140 // For each offset, look up the value in the map and only evaluate function if it's not there
141 for (int offset=1; offset<5; offset++) {
142 int location = n-offset;
143 std::map<int, std::complex<Real> >::const_iterator position = fLookupTable_.find(location);
144 if (position != fLookupTable_.end()) {
145 std::complex<Real> value = position->second;
146 fMinusN[offset-1] = value;
147 } else {
148 fMinusN[offset-1] = f(z1, z2, z3, z4, location, tau);
149 }
150 }
151
152 result = prefactor * (z1*tau*tau*fMinusN[3] + z2*tau*fMinusN[2] + (z3 - 0.5*kappa_*kappa_/(sigma_*sigma_))*fMinusN[1]);
153 }
154
155 // Store this value in our mutable lookup map
156 fLookupTable_[n] = result;
157
158 return result;
159 };
160
161 std::pair<std::complex<Real>, std::complex<Real> >
163 const std::complex<Real>& z1,
164 const std::complex<Real>& z2,
165 const std::complex<Real>& z3,
166 const std::complex<Real>& z4,
167 Real tau,
168 Size cutoff) const {
169 std::complex<Real> temp = 0.0;
170 std::complex<Real> runningSum1 = 0.0;
171 std::complex<Real> runningSum2 = 0.0;
172
173 for (Size i=0; i<cutoff; i++) {
174 temp = f(z1, z2, z3, z4, i, tau);
175 runningSum1 += temp;
176 runningSum2 += temp*Real(i)/tau;
177 }
178
179 std::pair<std::complex<Real>, std::complex<Real> > result(runningSum1, runningSum2);
180
181 return result;
182 };
183
185 const std::complex<Real>& s,
186 const std::complex<Real>& w,
187 Real T,
188 Real t,
189 Size cutoff) const {
190 Real tau = T - t;
191
192 std::complex<Real> z1 = z1_f(s, w, T);
193 std::complex<Real> z2 = z2_f(s, w, T);
194 std::complex<Real> z3 = z3_f(s, w, T);
195 std::complex<Real> z4 = z4_f(s, w);
196
197 // Clear the mutable lookup map before calling fLookupTable
198 fLookupTable_ = std::map<int, std::complex<Real> >();
199 std::pair<std::complex<Real>, std::complex<Real> > temp = F_F_tilde(z1, z2, z3, z4, tau, cutoff);
200
201 std::complex<Real> F, F_tilde;
202 F = temp.first;
203 F_tilde = temp.second;
204
205 return std::exp(-a1_*F_tilde/F - a2_*std::log(F) + a3_*s + a4_*w + a5_);
206 }
207
210 "not a geometric average option");
211 QL_REQUIRE(arguments_.exercise->type() == Exercise::European,
212 "not an European Option");
213
214 ext::shared_ptr<PlainVanillaPayoff> payoff =
215 ext::dynamic_pointer_cast<PlainVanillaPayoff>(arguments_.payoff);
216 QL_REQUIRE(payoff, "non-plain payoff given");
217
218 Real strike = payoff->strike();
219 Date exercise = arguments_.exercise->lastDate();
220
221 Time expiryTime = this->process_->time(exercise);
222 QL_REQUIRE(expiryTime >= 0.0, "Expiry Date cannot be in the past");
223
224 Real expiryDcf = riskFreeRate_->discount(expiryTime);
225 Real expiryDividendDiscount = dividendYield_->discount(expiryTime);
226
227 // TODO: extend to cover seasoned options (discussed in paper)
228 Time startTime = 0.0;
229
230 // These parameters only need to be calculated once per pricing, but are
231 // functions of t and T so need to be reset in calculate()
232 Time t = startTime;
233 Time T = expiryTime;
234 Time tau = T - t;
235 Real logS0 = std::log(s0_->value());
236
237 // To deal with non-constant rates and dividends, we reformulate Eq.s (14) to (17) with
238 // r_ --> (r(t) - q(t)), which gives the new expressions for a3 and a4 used below
239 Real dcf = riskFreeRate_->discount(T) / riskFreeRate_->discount(t);
240 Real qdcf = dividendYield_->discount(T) / dividendYield_->discount(t);
241 DcfIntegrand dcfIntegrand = DcfIntegrand(t, T, riskFreeRate_, dividendYield_);
242 Real integratedDcf = integrator_(dcfIntegrand);
243
244 a3_ = (tau*logS0 + integratedDcf)/T - kappa_*theta_*rho_*tau*tau/(2*sigma_*T) - rho_*tau*v0_/(sigma_*T);
245 a4_ = logS0*qdcf/dcf - rho_*v0_/sigma_ + rho_*kappa_*theta_*tau/sigma_;
247
248 // Calculate the two terms in eq (29) - Phi(1,0) is real (asian forward) but need to type convert
249 Real term1 = 0.5 * (std::real(Phi(1,0, T, t, summationCutoff_)) - strike);
250
251 Integrand integrand(T, summationCutoff_, strike, this, xiRightLimit_);
252 Real term2 = integrator_(integrand) / M_PI;
253
254 // Apply the payoff functions
255 Real value = 0.0;
256 switch (payoff->optionType()){
257 case Option::Call:
258 value = expiryDcf * (term1 + term2);
259 break;
260 case Option::Put:
261 value = expiryDcf * (-term1 + term2);
262 break;
263 default:
264 QL_FAIL("unknown option type");
265 }
266
267 results_.value = value;
268
269 results_.additionalResults["dcf"] = expiryDcf;
270 results_.additionalResults["qf"] = expiryDividendDiscount;
271 results_.additionalResults["s0"] = s0_->value();
272 results_.additionalResults["strike"] = strike;
273 results_.additionalResults["expiryTime"] = expiryTime;
274 results_.additionalResults["exercise"] = exercise;
275
276 results_.additionalResults["term1"] = term1;
277 results_.additionalResults["term2"] = term2;
278 results_.additionalResults["xiRightLimit"] = xiRightLimit_;
279 results_.additionalResults["summationCutoff"] = summationCutoff_;
280
286 }
287}
Analytic engine for continuous geometric average price Asian in the Heston model.
std::complex< Real > f(const std::complex< Real > &z1, const std::complex< Real > &z2, const std::complex< Real > &z3, const std::complex< Real > &z4, int n, Real tau) const
std::complex< Real > z3_f(const std::complex< Real > &s, const std::complex< Real > &w, Real T) const
std::complex< Real > z4_f(const std::complex< Real > &s, const std::complex< Real > &w) const
std::pair< std::complex< Real >, std::complex< Real > > F_F_tilde(const std::complex< Real > &z1, const std::complex< Real > &z2, const std::complex< Real > &z3, const std::complex< Real > &z4, Real tau, Size cutoff=50) const
std::complex< Real > z2_f(const std::complex< Real > &s, const std::complex< Real > &w, Real T) const
std::complex< Real > z1_f(const std::complex< Real > &s, const std::complex< Real > &w, Real T) const
std::complex< Real > Phi(const std::complex< Real > &s, const std::complex< Real > &w, Real T, Real t=0.0, Size cutoff=50) const
AnalyticContinuousGeometricAveragePriceAsianHestonEngine(ext::shared_ptr< HestonProcess > process, Size summationCutoff=50, Real xiRightLimit=100.0)
Concrete date class.
Definition: date.hpp:125
std::map< std::string, ext::any > additionalResults
Definition: instrument.hpp:123
std::pair< iterator, bool > registerWith(const ext::shared_ptr< Observable > &)
Definition: observable.hpp:228
const DefaultType & t
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
Definition: errors.hpp:117
#define QL_FAIL(message)
throw an error (possibly with file and line information)
Definition: errors.hpp:92
Real Time
continuous quantity with 1-year units
Definition: types.hpp:62
QL_REAL Real
real number
Definition: types.hpp:50
std::size_t Size
size of a container
Definition: types.hpp:58
ext::shared_ptr< QuantLib::Payoff > payoff
#define M_PI
Definition: any.hpp:35
STL namespace.
Real F
Definition: sabr.cpp:200