QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
gsrprocesscore.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) 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
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
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(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
69 const Time dt) const {
70 Real t = w + dt;
71 std::pair<Real, Real> key;
72 key = std::make_pair(w, t);
73 std::map<std::pair<Real, Real>, Real>::const_iterator k = cache1_.find(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(w); i <= upperIndex(t) - 1; i++) {
79 res2 *= exp(-rev(i) * (cappedTime(i + 1, t) - flooredTime(i, w)));
80 }
81 cache1_.insert(std::make_pair(key, res2));
82 return res2 * xw;
83}
84
86 const Time dt) const {
87
88 Real t = w + dt;
89
90 std::pair<Real, Real> key;
91 key = std::make_pair(w, t);
92 std::map<std::pair<Real, Real>, Real>::const_iterator k =
93 cache2a_.find(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(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(l) ? Real(vol(l) * vol(l) * (time2(l + 1) - time2(l)))
106 : vol(l) * vol(l) / (2.0 * rev(l)) *
107 (1.0 - exp(-2.0 * rev(l) *
108 (time2(l + 1) - time2(l))));
109 // zeta_i (i>k)
110 for (int i = k + 1; i <= upperIndex(t) - 1; i++)
111 res2 *= exp(-rev(i) * (cappedTime(i + 1, t) - time2(i)));
112 // beta_j (j<k)
113 for (int j = l + 1; j <= k - 1; j++)
114 res2 *= exp(-2.0 * rev(j) * (time2(j + 1) - time2(j)));
115 // zeta_k beta_k
116 res2 *=
117 revZero(k)
118 ? Real(2.0 * time2(k) - flooredTime(k, w) -
119 cappedTime(k + 1, t) -
120 2.0 * (time2(k) - cappedTime(k + 1, t)))
121 : Real((exp(rev(k) * (2.0 * time2(k) - flooredTime(k, w) -
122 cappedTime(k + 1, t))) -
123 exp(2.0 * rev(k) * (time2(k) - cappedTime(k + 1, t)))) /
124 rev(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(k)
133 ? Real(vol(k) * vol(k) / 4.0 *
134 (4.0 * pow(cappedTime(k + 1, t) - time2(k), 2.0) -
135 (pow(flooredTime(k, w) - 2.0 * time2(k) +
136 cappedTime(k + 1, t),
137 2.0) +
138 pow(cappedTime(k + 1, t) - flooredTime(k, w), 2.0))))
139 : Real(vol(k) * vol(k) / (2.0 * rev(k) * rev(k)) *
140 (exp(-2.0 * rev(k) * (cappedTime(k + 1, t) - time2(k))) +
141 1.0 -
142 (exp(-rev(k) * (flooredTime(k, w) - 2.0 * time2(k) +
143 cappedTime(k + 1, t))) +
144 exp(-rev(k) *
145 (cappedTime(k + 1, t) - flooredTime(k, w))))));
146 // zeta_i (i>k)
147 for (int i = k + 1; i <= upperIndex(t) - 1; i++)
148 res2 *= exp(-rev(i) * (cappedTime(i + 1, t) - time2(i)));
149 // no beta_j in this case ...
150 res += res2;
151 }
152
153 cache2a_.insert(std::make_pair(key, res));
154
155 return res;
156} // expectation_rn_part
157
159 const Time dt) const {
160
161 Real t = w + dt;
162
163 std::pair<Real, Real> key;
164 key = std::make_pair(w, t);
165 std::map<std::pair<Real, Real>, Real>::const_iterator k =
166 cache2b_.find(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(w); k <= upperIndex(t) - 1; k++) {
173 Real res2 = 0.0;
174 // l>k
175 for (int l = k + 1; l <= upperIndex(T_) - 1; l++) {
176 Real res3 = 1.0;
177 // eta_l
178 res3 *= revZero(l)
179 ? Real(cappedTime(l + 1, T_) - time2(l))
180 : (1.0 -
181 exp(-rev(l) * (cappedTime(l + 1, T_) - time2(l)))) /
182 rev(l);
183 // zeta_i (i>k)
184 for (int i = k + 1; i <= upperIndex(t) - 1; i++)
185 res3 *= exp(-rev(i) * (cappedTime(i + 1, t) - time2(i)));
186 // gamma_j (j>k)
187 for (int j = k + 1; j <= l - 1; j++)
188 res3 *= exp(-rev(j) * (time2(j + 1) - time2(j)));
189 // zeta_k gamma_k
190 res3 *=
191 revZero(k)
192 ? Real((cappedTime(k + 1, t) - time2(k + 1) -
193 (2.0 * flooredTime(k, w) - cappedTime(k + 1, t) -
194 time2(k + 1))) /
195 2.0)
196 : Real((exp(rev(k) * (cappedTime(k + 1, t) - time2(k + 1))) -
197 exp(rev(k) * (2.0 * flooredTime(k, w) -
198 cappedTime(k + 1, t) - time2(k + 1)))) /
199 (2.0 * rev(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(k)
208 ? Real((-pow(cappedTime(k + 1, t) - cappedTime(k + 1, T_), 2.0) -
209 2.0 * pow(cappedTime(k + 1, t) - flooredTime(k, w), 2.0) +
210 pow(2.0 * flooredTime(k, w) - cappedTime(k + 1, T_) -
211 cappedTime(k + 1, t),
212 2.0)) /
213 4.0)
214 : Real((2.0 - exp(rev(k) *
215 (cappedTime(k + 1, t) - cappedTime(k + 1, T_))) -
216 (2.0 * exp(-rev(k) *
217 (cappedTime(k + 1, t) - flooredTime(k, w))) -
218 exp(rev(k) *
219 (2.0 * flooredTime(k, w) - cappedTime(k + 1, T_) -
220 cappedTime(k + 1, t))))) /
221 (2.0 * rev(k) * rev(k)));
222 // zeta_i (i>k)
223 for (int i = k + 1; i <= upperIndex(t) - 1; i++)
224 res3 *= exp(-rev(i) * (cappedTime(i + 1, t) - time2(i)));
225 // no gamma_j in this case ...
226 res2 += res3;
227 // add to main accumulator
228 res += -vol(k) * vol(k) * res2;
229 }
230
231 cache2b_.insert(std::make_pair(key, 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(w, t);
242 std::map<std::pair<Real, Real>, Real>::const_iterator k = cache3_.find(key);
243 if (k != cache3_.end())
244 return k->second;
245
246 Real res = 0.0;
247 for (int k = lowerIndex(w); k <= upperIndex(t) - 1; k++) {
248 Real res2 = vol(k) * vol(k);
249 // zeta_k^2
250 res2 *= revZero(k)
251 ? Real(-(flooredTime(k, w) - cappedTime(k + 1, t)))
252 : (1.0 - exp(2.0 * rev(k) *
253 (flooredTime(k, w) - cappedTime(k + 1, t)))) /
254 (2.0 * rev(k));
255 // zeta_i (i>k)
256 for (int i = k + 1; i <= upperIndex(t) - 1; i++) {
257 res2 *= exp(-2.0 * rev(i) * (cappedTime(i + 1, t) - time2(i)));
258 }
259 res += res2;
260 }
261
262 cache3_.insert(std::make_pair(key, res));
263 return res;
264}
265
267 Real key;
268 key = t;
269 std::map<Real, Real>::const_iterator k = cache4_.find(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(-2.0 * rev(j) * (cappedTime(j + 1, t) - time2(j)));
278 }
279 res2 *= revZero(i) ? Real(vol(i) * vol(i) * (cappedTime(i + 1, t) - time2(i)))
280 : (vol(i) * vol(i) / (2.0 * rev(i)) *
281 (1.0 - exp(-2.0 * rev(i) *
282 (cappedTime(i + 1, t) - time2(i)))));
283 res += res2;
284 }
285
286 cache4_.insert(std::make_pair(key, 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(w, t);
293 std::map<std::pair<Real, Real>, Real>::const_iterator k = cache5_.find(key);
294 if (k != cache5_.end())
295 return k->second;
296
297 Real res = 0.0;
298 for (int i = lowerIndex(t); i <= upperIndex(w) - 1; i++) {
299 Real res2 = 1.0;
300 for (int j = lowerIndex(t); j <= i - 1; j++) {
301 res2 *= exp(-rev(j) * (time2(j + 1) - flooredTime(j, t)));
302 }
303 res2 *= revZero(i) ? Real(cappedTime(i + 1, w) - flooredTime(i, t))
304 : (1.0 - exp(-rev(i) * (cappedTime(i + 1, w) -
305 flooredTime(i, t)))) /
306 rev(i);
307 res += res2;
308 }
309
310 cache5_.insert(std::make_pair(key, res));
311 return res;
312}
313
315 return static_cast<int>(std::upper_bound(times_.begin(), times_.end(), t) -
316 times_.begin());
317}
318
321 return 0;
322 return static_cast<int>(
323 std::upper_bound(times_.begin(), times_.end(), 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(cap, time2(index)) : time2(index);
330}
331
333 const Real floor) const {
334 return floor != Null<Real>() ? std::max(floor, 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
1-D array used in linear algebra.
Definition: array.hpp:52
const_iterator end() const
Definition: array.hpp:511
Real back() const
Definition: array.hpp:458
Size size() const
dimension of the array
Definition: array.hpp:495
const_iterator begin() const
Definition: array.hpp:503
template class providing a null value for a given type.
Definition: null.hpp:76
Real cappedTime(Size index, Real cap=Null< Real >()) const
GsrProcessCore(const Array &times, const Array &vols, const Array &reversions, Real T=60.0)
Real expectation_x0dep_part(Time w, Real xw, Time dt) const
std::map< Real, Real > cache4_
Real expectation_rn_part(Time w, Time dt) const
Real flooredTime(Size index, Real floor=Null< Real >()) const
std::map< std::pair< Real, Real >, Real > cache2a_
Real expectation_tf_part(Time w, Time dt) const
std::map< std::pair< Real, Real >, Real > cache1_
bool revZero(Size index) const
std::map< std::pair< Real, Real >, Real > cache2b_
std::map< std::pair< Real, Real >, Real > cache3_
Real G(Time t, Time w) const
std::map< std::pair< Real, Real >, Real > cache5_
Real variance(Time w, Time dt) const
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_EPSILON
Definition: qldefines.hpp:178
#define QL_MIN_POSITIVE_REAL
Definition: qldefines.hpp:177
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
Core computations for the gsr process in risk neutral and T-forward measure.
Definition: any.hpp:35