1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2015 Peter Caspers
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/processes/gsrprocesscore.hpp>
21#include <cmath>
22
23using std::exp;
24using std::pow;
25
26namespace QuantLib {
27
28namespace detail {
29
30GsrProcessCore::GsrProcessCore(const Array &times, const Array &vols,
31 const Array &reversions, const Real T)
32 : times_(times), vols_(vols), reversions_(reversions),
33 T_(T), revZero_(reversions.size(), false) {
34
35 QL_REQUIRE(times.size() == vols.size() - 1,
36 "number of volatilities ("
37 << vols.size() << ") compared to number of times ("
38 << times_.size() << " must be bigger by one");
39 QL_REQUIRE(times.size() == reversions.size() - 1 || reversions.size() == 1,
40 "number of reversions ("
41 << vols.size() << ") compared to number of times ("
42 << times_.size() << " must be bigger by one, or exactly "
43 "1 reversion must be given");
44 for (int i = 0; i < ((int)times.size()) - 1; i++)
45 QL_REQUIRE(times[i] < times[i + 1], "times must be increasing ("
46 << times[i] << "@" << i << " , "
47 << times[i + 1] << "@" << i + 1
48 << ")");
49 flushCache();
50}
51
52void GsrProcessCore::flushCache() const {
53 for (int i = 0; i < (int)reversions_.size(); i++)
54 // small reversions cause numerical problems, so we keep them
55 // away from zero
56 if (std::fabs(x: reversions_[i]) < 1E-4)
57 revZero_[i] = true;
58 else
59 revZero_[i] = false;
60 cache1_.clear();
61 cache2a_.clear();
62 cache2b_.clear();
63 cache3_.clear();
64 cache4_.clear();
65 cache5_.clear();
66}
67
68Real GsrProcessCore::expectation_x0dep_part(const Time w, const Real xw,
69 const Time dt) const {
70 Real t = w + dt;
71 std::pair<Real, Real> key;
72 key = std::make_pair(x: w, y&: t);
73 std::map<std::pair<Real, Real>, Real>::const_iterator k = cache1_.find(x: key);
74 if (k != cache1_.end())
75 return xw * (k->second);
76 // A(w,t)x(w)
77 Real res2 = 1.0;
78 for (int i = lowerIndex(t: w); i <= upperIndex(t) - 1; i++) {
79 res2 *= exp(x: -rev(index: i) * (cappedTime(index: i + 1, cap: t) - flooredTime(index: i, floor: w)));
80 }
81 cache1_.insert(x: std::make_pair(x&: key, y&: res2));
82 return res2 * xw;
83}
84
85Real GsrProcessCore::expectation_rn_part(const Time w,
86 const Time dt) const {
87
88 Real t = w + dt;
89
90 std::pair<Real, Real> key;
91 key = std::make_pair(x: w, y&: t);
92 std::map<std::pair<Real, Real>, Real>::const_iterator k =
93 cache2a_.find(x: key);
94 if (k != cache2a_.end())
95 return k->second;
96
97 Real res = 0.0;
98
99 // \int A(s,t)y(s)
100 for (int k = lowerIndex(t: w); k <= upperIndex(t) - 1; k++) {
101 // l<k
102 for (int l = 0; l <= k - 1; l++) {
103 Real res2 = 1.0;
104 // alpha_l
105 res2 *= revZero(index: l) ? Real(vol(index: l) * vol(index: l) * (time2(index: l + 1) - time2(index: l)))
106 : vol(index: l) * vol(index: l) / (2.0 * rev(index: l)) *
107 (1.0 - exp(x: -2.0 * rev(index: l) *
108 (time2(index: l + 1) - time2(index: l))));
109 // zeta_i (i>k)
110 for (int i = k + 1; i <= upperIndex(t) - 1; i++)
111 res2 *= exp(x: -rev(index: i) * (cappedTime(index: i + 1, cap: t) - time2(index: i)));
112 // beta_j (j<k)
113 for (int j = l + 1; j <= k - 1; j++)
114 res2 *= exp(x: -2.0 * rev(index: j) * (time2(index: j + 1) - time2(index: j)));
115 // zeta_k beta_k
116 res2 *=
117 revZero(index: k)
118 ? Real(2.0 * time2(index: k) - flooredTime(index: k, floor: w) -
119 cappedTime(index: k + 1, cap: t) -
120 2.0 * (time2(index: k) - cappedTime(index: k + 1, cap: t)))
121 : Real((exp(x: rev(index: k) * (2.0 * time2(index: k) - flooredTime(index: k, floor: w) -
122 cappedTime(index: k + 1, cap: t))) -
123 exp(x: 2.0 * rev(index: k) * (time2(index: k) - cappedTime(index: k + 1, cap: t)))) /
124 rev(index: k));
125 // add to sum
126 res += res2;
127 }
128 // l=k
129 Real res2 = 1.0;
130 // alpha_k zeta_k
131 res2 *=
132 revZero(index: k)
133 ? Real(vol(index: k) * vol(index: k) / 4.0 *
134 (4.0 * pow(x: cappedTime(index: k + 1, cap: t) - time2(index: k), y: 2.0) -
135 (pow(x: flooredTime(index: k, floor: w) - 2.0 * time2(index: k) +
136 cappedTime(index: k + 1, cap: t),
137 y: 2.0) +
138 pow(x: cappedTime(index: k + 1, cap: t) - flooredTime(index: k, floor: w), y: 2.0))))
139 : Real(vol(index: k) * vol(index: k) / (2.0 * rev(index: k) * rev(index: k)) *
140 (exp(x: -2.0 * rev(index: k) * (cappedTime(index: k + 1, cap: t) - time2(index: k))) +
141 1.0 -
142 (exp(x: -rev(index: k) * (flooredTime(index: k, floor: w) - 2.0 * time2(index: k) +
143 cappedTime(index: k + 1, cap: t))) +
144 exp(x: -rev(index: k) *
145 (cappedTime(index: k + 1, cap: t) - flooredTime(index: k, floor: w))))));
146 // zeta_i (i>k)
147 for (int i = k + 1; i <= upperIndex(t) - 1; i++)
148 res2 *= exp(x: -rev(index: i) * (cappedTime(index: i + 1, cap: t) - time2(index: i)));
149 // no beta_j in this case ...
150 res += res2;
151 }
152
153 cache2a_.insert(x: std::make_pair(x&: key, y&: res));
154
155 return res;
156} // expectation_rn_part
157
158Real GsrProcessCore::expectation_tf_part(const Time w,
159 const Time dt) const {
160
161 Real t = w + dt;
162
163 std::pair<Real, Real> key;
164 key = std::make_pair(x: w, y&: t);
165 std::map<std::pair<Real, Real>, Real>::const_iterator k =
166 cache2b_.find(x: key);
167 if (k != cache2b_.end())
168 return k->second;
169
170 Real res = 0.0;
171 // int -A(s,t) \sigma^2 G(s,T)
172 for (int k = lowerIndex(t: w); k <= upperIndex(t) - 1; k++) {
173 Real res2 = 0.0;
174 // l>k
175 for (int l = k + 1; l <= upperIndex(t: T_) - 1; l++) {
176 Real res3 = 1.0;
177 // eta_l
178 res3 *= revZero(index: l)
179 ? Real(cappedTime(index: l + 1, cap: T_) - time2(index: l))
180 : (1.0 -
181 exp(x: -rev(index: l) * (cappedTime(index: l + 1, cap: T_) - time2(index: l)))) /
182 rev(index: l);
183 // zeta_i (i>k)
184 for (int i = k + 1; i <= upperIndex(t) - 1; i++)
185 res3 *= exp(x: -rev(index: i) * (cappedTime(index: i + 1, cap: t) - time2(index: i)));
186 // gamma_j (j>k)
187 for (int j = k + 1; j <= l - 1; j++)
188 res3 *= exp(x: -rev(index: j) * (time2(index: j + 1) - time2(index: j)));
189 // zeta_k gamma_k
190 res3 *=
191 revZero(index: k)
192 ? Real((cappedTime(index: k + 1, cap: t) - time2(index: k + 1) -
193 (2.0 * flooredTime(index: k, floor: w) - cappedTime(index: k + 1, cap: t) -
194 time2(index: k + 1))) /
195 2.0)
196 : Real((exp(x: rev(index: k) * (cappedTime(index: k + 1, cap: t) - time2(index: k + 1))) -
197 exp(x: rev(index: k) * (2.0 * flooredTime(index: k, floor: w) -
198 cappedTime(index: k + 1, cap: t) - time2(index: k + 1)))) /
199 (2.0 * rev(index: k)));
200 // add to sum
201 res2 += res3;
202 }
203 // l=k
204 Real res3 = 1.0;
205 // eta_k zeta_k
206 res3 *=
207 revZero(index: k)
208 ? Real((-pow(x: cappedTime(index: k + 1, cap: t) - cappedTime(index: k + 1, cap: T_), y: 2.0) -
209 2.0 * pow(x: cappedTime(index: k + 1, cap: t) - flooredTime(index: k, floor: w), y: 2.0) +
210 pow(x: 2.0 * flooredTime(index: k, floor: w) - cappedTime(index: k + 1, cap: T_) -
211 cappedTime(index: k + 1, cap: t),
212 y: 2.0)) /
213 4.0)
214 : Real((2.0 - exp(x: rev(index: k) *
215 (cappedTime(index: k + 1, cap: t) - cappedTime(index: k + 1, cap: T_))) -
216 (2.0 * exp(x: -rev(index: k) *
217 (cappedTime(index: k + 1, cap: t) - flooredTime(index: k, floor: w))) -
218 exp(x: rev(index: k) *
219 (2.0 * flooredTime(index: k, floor: w) - cappedTime(index: k + 1, cap: T_) -
220 cappedTime(index: k + 1, cap: t))))) /
221 (2.0 * rev(index: k) * rev(index: k)));
222 // zeta_i (i>k)
223 for (int i = k + 1; i <= upperIndex(t) - 1; i++)
224 res3 *= exp(x: -rev(index: i) * (cappedTime(index: i + 1, cap: t) - time2(index: i)));
225 // no gamma_j in this case ...
226 res2 += res3;
227 // add to main accumulator
228 res += -vol(index: k) * vol(index: k) * res2;
229 }
230
231 cache2b_.insert(x: std::make_pair(x&: key, y&: res));
232
233 return res;
234} // expectation_tf_part
235
236Real GsrProcessCore::variance(const Time w, const Time dt) const {
237
238 Real t = w + dt;
239
240 std::pair<Real, Real> key;
241 key = std::make_pair(x: w, y&: t);
242 std::map<std::pair<Real, Real>, Real>::const_iterator k = cache3_.find(x: key);
243 if (k != cache3_.end())
244 return k->second;
245
246 Real res = 0.0;
247 for (int k = lowerIndex(t: w); k <= upperIndex(t) - 1; k++) {
248 Real res2 = vol(index: k) * vol(index: k);
249 // zeta_k^2
250 res2 *= revZero(index: k)
251 ? Real(-(flooredTime(index: k, floor: w) - cappedTime(index: k + 1, cap: t)))
252 : (1.0 - exp(x: 2.0 * rev(index: k) *
253 (flooredTime(index: k, floor: w) - cappedTime(index: k + 1, cap: t)))) /
254 (2.0 * rev(index: k));
255 // zeta_i (i>k)
256 for (int i = k + 1; i <= upperIndex(t) - 1; i++) {
257 res2 *= exp(x: -2.0 * rev(index: i) * (cappedTime(index: i + 1, cap: t) - time2(index: i)));
258 }
259 res += res2;
260 }
261
262 cache3_.insert(x: std::make_pair(x&: key, y&: res));
263 return res;
264}
265
266Real GsrProcessCore::y(const Time t) const {
267 Real key;
268 key = t;
269 std::map<Real, Real>::const_iterator k = cache4_.find(x: key);
270 if (k != cache4_.end())
271 return k->second;
272
273 Real res = 0.0;
274 for (int i = 0; i <= upperIndex(t) - 1; i++) {
275 Real res2 = 1.0;
276 for (int j = i + 1; j <= upperIndex(t) - 1; j++) {
277 res2 *= exp(x: -2.0 * rev(index: j) * (cappedTime(index: j + 1, cap: t) - time2(index: j)));
278 }
279 res2 *= revZero(index: i) ? Real(vol(index: i) * vol(index: i) * (cappedTime(index: i + 1, cap: t) - time2(index: i)))
280 : (vol(index: i) * vol(index: i) / (2.0 * rev(index: i)) *
281 (1.0 - exp(x: -2.0 * rev(index: i) *
282 (cappedTime(index: i + 1, cap: t) - time2(index: i)))));
283 res += res2;
284 }
285
286 cache4_.insert(x: std::make_pair(x&: key, y&: res));
287 return res;
288}
289
290Real GsrProcessCore::G(const Time t, const Time w) const {
291 std::pair<Real, Real> key;
292 key = std::make_pair(x: w, y: t);
293 std::map<std::pair<Real, Real>, Real>::const_iterator k = cache5_.find(x: key);
294 if (k != cache5_.end())
295 return k->second;
296
297 Real res = 0.0;
298 for (int i = lowerIndex(t); i <= upperIndex(t: w) - 1; i++) {
299 Real res2 = 1.0;
300 for (int j = lowerIndex(t); j <= i - 1; j++) {
301 res2 *= exp(x: -rev(index: j) * (time2(index: j + 1) - flooredTime(index: j, floor: t)));
302 }
303 res2 *= revZero(index: i) ? Real(cappedTime(index: i + 1, cap: w) - flooredTime(index: i, floor: t))
304 : (1.0 - exp(x: -rev(index: i) * (cappedTime(index: i + 1, cap: w) -
305 flooredTime(index: i, floor: t)))) /
306 rev(index: i);
307 res += res2;
308 }
309
310 cache5_.insert(x: std::make_pair(x&: key, y&: res));
311 return res;
312}
313
314int GsrProcessCore::lowerIndex(const Time t) const {
315 return static_cast<int>(std::upper_bound(first: times_.begin(), last: times_.end(), val: t) -
316 times_.begin());
317}
318
319int GsrProcessCore::upperIndex(const Time t) const {
320 if (t < QL_MIN_POSITIVE_REAL)
321 return 0;
322 return static_cast<int>(
323 std::upper_bound(first: times_.begin(), last: times_.end(), val: t - QL_EPSILON) -
324 times_.begin()) +
325 1;
326}
327
328Real GsrProcessCore::cappedTime(const Size index, const Real cap) const {
329 return cap != Null<Real>() ? std::min(a: cap, b: time2(index)) : time2(index);
330}
331
332Real GsrProcessCore::flooredTime(const Size index,
333 const Real floor) const {
334 return floor != Null<Real>() ? std::max(a: floor, b: time2(index)) : time2(index);
335}
336
337Real GsrProcessCore::time2(const Size index) const {
338 if (index == 0)
339 return 0.0;
340 if (index > times_.size())
341 return T_; // FIXME how to ensure that forward
342 // measure time is geq all times
343 // given
344 return times_[index - 1];
345}
346
347Real GsrProcessCore::vol(const Size index) const {
348 if (index >= vols_.size())
349 return vols_.back();
350 return vols_[index];
351}
352
353Real GsrProcessCore::rev(const Size index) const {
354 if (index >= reversions_.size())
355 return reversions_.back();
356 return reversions_[index];
357}
358
359bool GsrProcessCore::revZero(const Size index) const {
360 if (index >= revZero_.size())
361 return revZero_.back();
362 return revZero_[index];
363}
364
365} // namespace detail
366
367} // namesapce QuantLib
368

source code of quantlib/ql/processes/gsrprocesscore.cpp