1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2001, 2002, 2003 Sadruddin Rejeb
5 Copyright (C) 2004 Mike Parker
6 Copyright (C) 2021 Magnus Mencke
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 <ql/math/distributions/normaldistribution.hpp>
23#include <ql/math/integrals/segmentintegral.hpp>
24#include <ql/math/solvers1d/brent.hpp>
25#include <ql/models/shortrate/twofactormodels/g2.hpp>
26#include <ql/pricingengines/blackformula.hpp>
27#include <utility>
28
29namespace QuantLib {
30
31 G2::G2(const Handle<YieldTermStructure>& termStructure,
32 Real a, Real sigma, Real b, Real eta, Real rho)
33 : TwoFactorModel(5), TermStructureConsistentModel(termStructure),
34 a_(arguments_[0]), sigma_(arguments_[1]), b_(arguments_[2]),
35 eta_(arguments_[3]), rho_(arguments_[4]) {
36
37 a_ = ConstantParameter(a, PositiveConstraint());
38 sigma_ = ConstantParameter(sigma, PositiveConstraint());
39 b_ = ConstantParameter(b, PositiveConstraint());
40 eta_ = ConstantParameter(eta, PositiveConstraint());
41 rho_ = ConstantParameter(rho, BoundaryConstraint(-1.0, 1.0));
42
43 G2::generateArguments();
44
45 registerWith(h: termStructure);
46 }
47
48 ext::shared_ptr<TwoFactorModel::ShortRateDynamics> G2::dynamics() const {
49 return ext::shared_ptr<ShortRateDynamics>(new
50 Dynamics(phi_, a(), sigma(), b(), eta(), rho()));
51 }
52
53 void G2::generateArguments() {
54
55 phi_ = FittingParameter(termStructure(),
56 a(), sigma(), b(), eta(), rho());
57 }
58
59 Real G2::sigmaP(Time t, Time s) const {
60 Real temp = 1.0 - std::exp(x: -(a()+b())*t);
61 Real temp1 = 1.0 - std::exp(x: -a()*(s-t));
62 Real temp2 = 1.0 - std::exp(x: -b()*(s-t));
63 Real a3 = a()*a()*a();
64 Real b3 = b()*b()*b();
65 Real sigma2 = sigma()*sigma();
66 Real eta2 = eta()*eta();
67 Real value =
68 0.5*sigma2*temp1*temp1*(1.0 - std::exp(x: -2.0*a()*t))/a3 +
69 0.5*eta2*temp2*temp2*(1.0 - std::exp(x: -2.0*b()*t))/b3 +
70 2.0*rho()*sigma()*eta()/(a()*b()*(a()+b()))*
71 temp1*temp2*temp;
72 return std::sqrt(x: value);
73 }
74
75 Real G2::discountBond(Time t, Time T, Real x, Real y) const {
76 return A(t,T) * std::exp(x: -B(x: a(),t: (T-t))*x-B(x: b(),t: (T-t))*y);
77 }
78
79 Real G2::discountBondOption(Option::Type type, Real strike, Time maturity,
80 Time bondMaturity) const {
81
82 Real v = sigmaP(t: maturity, s: bondMaturity);
83 Real f = termStructure()->discount(t: bondMaturity);
84 Real k = termStructure()->discount(t: maturity)*strike;
85
86 return blackFormula(optionType: type, strike: k, forward: f, stdDev: v);
87 }
88
89 Real G2::V(Time t) const {
90 Real expat = std::exp(x: -a()*t);
91 Real expbt = std::exp(x: -b()*t);
92 Real cx = sigma()/a();
93 Real cy = eta()/b();
94 Real valuex = cx*cx*(t + (2.0*expat-0.5*expat*expat-1.5)/a());
95 Real valuey = cy*cy*(t + (2.0*expbt-0.5*expbt*expbt-1.5)/b());
96 Real value = 2.0*rho()*cx*cy* (t + (expat - 1.0)/a()
97 + (expbt - 1.0)/b()
98 - (expat*expbt-1.0)/(a()+b()));
99 return valuex + valuey + value;
100 }
101
102 Real G2::A(Time t, Time T) const {
103 return termStructure()->discount(t: T)/termStructure()->discount(t)*
104 std::exp(x: 0.5*(V(t: T-t) - V(t: T) + V(t)));
105 }
106
107 Real G2::B(Real x, Time t) const {
108 return (1.0 - std::exp(x: -x*t))/x;
109 }
110
111 class G2::SwaptionPricingFunction {
112 public:
113 SwaptionPricingFunction(Real a,
114 Real sigma,
115 Real b,
116 Real eta,
117 Real rho,
118 Real w,
119 Real start,
120 std::vector<Time> payTimes,
121 Rate fixedRate,
122 const G2& model)
123 : a_(a), sigma_(sigma), b_(b), eta_(eta), rho_(rho), w_(w), T_(start),
124 t_(std::move(payTimes)), rate_(fixedRate), size_(t_.size()), A_(size_), Ba_(size_),
125 Bb_(size_) {
126
127
128 sigmax_ = sigma_*std::sqrt(x: 0.5*(1.0-std::exp(x: -2.0*a_*T_))/a_);
129 sigmay_ = eta_*std::sqrt(x: 0.5*(1.0-std::exp(x: -2.0*b_*T_))/b_);
130 rhoxy_ = rho_*eta_*sigma_*(1.0 - std::exp(x: -(a_+b_)*T_))/
131 ((a_+b_)*sigmax_*sigmay_);
132
133 Real temp = sigma_*sigma_/(a_*a_);
134 mux_ = -((temp+rho_*sigma_*eta_/(a_*b_))*(1.0 - std::exp(x: -a*T_)) -
135 0.5*temp*(1.0 - std::exp(x: -2.0*a_*T_)) -
136 rho_*sigma_*eta_/(b_*(a_+b_))*
137 (1.0- std::exp(x: -(b_+a_)*T_)));
138
139 temp = eta_*eta_/(b_*b_);
140 muy_ = -((temp+rho_*sigma_*eta_/(a_*b_))*(1.0 - std::exp(x: -b*T_)) -
141 0.5*temp*(1.0 - std::exp(x: -2.0*b_*T_)) -
142 rho_*sigma_*eta_/(a_*(a_+b_))*
143 (1.0- std::exp(x: -(b_+a_)*T_)));
144
145 for (Size i=0; i<size_; i++) {
146 A_[i] = model.A(t: T_, T: t_[i]);
147 Ba_[i] = model.B(x: a_, t: t_[i]-T_);
148 Bb_[i] = model.B(x: b_, t: t_[i]-T_);
149 }
150 }
151
152 Real mux() const { return mux_; }
153 Real sigmax() const { return sigmax_; }
154 Real operator()(Real x) const {
155 CumulativeNormalDistribution phi;
156 Real temp = (x - mux_)/sigmax_;
157 Real txy = std::sqrt(x: 1.0 - rhoxy_*rhoxy_);
158
159 Array lambda(size_);
160 Size i;
161 for (i=0; i<size_; i++) {
162 Real tau = (i==0 ? t_[0] - T_ : t_[i] - t_[i-1]);
163 Real c = (i==size_-1 ? Real(1.0+rate_*tau) : rate_*tau);
164 lambda[i] = c*A_[i]*std::exp(x: -Ba_[i]*x);
165 }
166
167 SolvingFunction function(lambda, Bb_) ;
168 Brent s1d;
169 s1d.setMaxEvaluations(1000);
170 Real searchBound = std::max(a: 10.0*sigmay_, b: 1.0);
171 Real yb = s1d.solve(f: function, accuracy: 1e-6, guess: 0.00, xMin: -searchBound, xMax: searchBound);
172
173 Real h1 = (yb - muy_)/(sigmay_*txy) -
174 rhoxy_*(x - mux_)/(sigmax_*txy);
175 Real value = phi(-w_*h1);
176
177
178 for (i=0; i<size_; i++) {
179 Real h2 = h1 +
180 Bb_[i]*sigmay_*std::sqrt(x: 1.0-rhoxy_*rhoxy_);
181 Real kappa = - Bb_[i] *
182 (muy_ - 0.5*txy*txy*sigmay_*sigmay_*Bb_[i] +
183 rhoxy_*sigmay_*(x-mux_)/sigmax_);
184 value -= lambda[i] *std::exp(x: kappa)*phi(-w_*h2);
185 }
186
187 return std::exp(x: -0.5*temp*temp)*value/
188 (sigmax_*std::sqrt(x: 2.0*M_PI));
189 }
190
191
192 private:
193 class SolvingFunction {
194 public:
195 SolvingFunction(const Array& lambda, const Array& Bb)
196 : lambda_(lambda), Bb_(Bb) {}
197 Real operator()(Real y) const {
198 Real value = 1.0;
199 for (Size i=0; i<lambda_.size(); i++) {
200 value -= lambda_[i]*std::exp(x: -Bb_[i]*y);
201 }
202 return value;
203 }
204 private:
205 const Array& lambda_;
206 const Array& Bb_;
207 };
208
209 Real a_, sigma_, b_, eta_, rho_, w_;
210 Time T_;
211 std::vector<Time> t_;
212 Rate rate_;
213 Size size_;
214 Array A_, Ba_, Bb_;
215 Real mux_, muy_, sigmax_, sigmay_, rhoxy_;
216 };
217
218 Real G2::swaption(const Swaption::arguments& arguments,
219 Rate fixedRate, Real range, Size intervals) const {
220
221 QL_REQUIRE(arguments.nominal != Null<Real>(),
222 "non-constant nominals are not supported yet");
223
224 Date settlement = termStructure()->referenceDate();
225 DayCounter dayCounter = termStructure()->dayCounter();
226 Time start = dayCounter.yearFraction(d1: settlement,
227 d2: arguments.floatingResetDates[0]);
228 Real w = (arguments.type==Swap::Payer ? 1 : -1 );
229
230 std::vector<Time> fixedPayTimes(arguments.fixedPayDates.size());
231 for (Size i=0; i<fixedPayTimes.size(); ++i)
232 fixedPayTimes[i] =
233 dayCounter.yearFraction(d1: settlement,
234 d2: arguments.fixedPayDates[i]);
235
236 SwaptionPricingFunction function(a(), sigma(), b(), eta(), rho(),
237 w, start,
238 fixedPayTimes,
239 fixedRate, (*this));
240
241 Real upper = function.mux() + range*function.sigmax();
242 Real lower = function.mux() - range*function.sigmax();
243 SegmentIntegral integrator(intervals);
244 return arguments.nominal * w * termStructure()->discount(t: start) *
245 integrator(function, lower, upper);
246 }
247
248}
249

source code of quantlib/ql/models/shortrate/twofactormodels/g2.cpp