QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
ratepseudorootjacobian.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
5Copyright (C) 2008 Mark Joshi
6
7This file is part of QuantLib, a free-software/open-source library
8for financial quantitative analysts and developers - http://quantlib.org/
9
10QuantLib is free software: you can redistribute it and/or modify it
11under the terms of the QuantLib license. You should have received a
12copy of the license along with this program; if not, please email
13<quantlib-dev@lists.sf.net>. The license is also available online at
14<http://quantlib.org/license.shtml>.
15
16This program is distributed in the hope that it will be useful, but WITHOUT
17ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18FOR A PARTICULAR PURPOSE. See the license for more details.
19*/
20
21
22
24#include <utility>
25
26namespace QuantLib
27{
28
29
31 Size aliveIndex,
32 Size numeraire,
33 const std::vector<Time>& taus,
34 const std::vector<Matrix>& pseudoBumps,
35 const std::vector<Spread>& displacements)
36 :
37 pseudoRoot_(pseudoRoot),
38 aliveIndex_(aliveIndex),
39 taus_(taus),
40 displacements_(displacements),
41 numberBumps_(pseudoBumps.size()),
42 factors_(pseudoRoot.columns()),
43 drifts_(taus.size()),
44 bumpedRates_(taus.size())
45 {
46 Size numberRates= taus.size();
47
48 QL_REQUIRE(pseudoRoot_.rows()==numberRates,
49 "pseudoRoot_.rows()<> taus.size()");
50
51 QL_REQUIRE(displacements_.size()==numberRates,
52 "displacements_.size()<> taus.size()");
53
54 QL_REQUIRE(drifts_.size()==numberRates,
55 "drifts_.size()<> taus.size()");
56
57 for (Size i=0; i < pseudoBumps.size(); ++i)
58 {
59 QL_REQUIRE(pseudoBumps[i].rows()==numberRates,
60 "pseudoBumps[i].rows()<> taus.size() with i =" << i);
61
62 QL_REQUIRE(pseudoBumps[i].columns()==factors_,
63 "pseudoBumps[i].columns()<> factors with i = " << i);
64
65
66 Matrix pseudo(pseudoRoot_);
67 pseudo += pseudoBumps[i];
68 pseudoBumped_.push_back(pseudo);
69 driftsComputers_.emplace_back(pseudo, displacements, taus, numeraire, aliveIndex);
70 }
71
72 }
73
74
75 void RatePseudoRootJacobianNumerical::getBumps(const std::vector<Rate>& oldRates,
76 const std::vector<Real>& , // not used in the numerical implementation
77 const std::vector<Rate>& newRates,
78 const std::vector<Real>& gaussians,
79 Matrix& B)
80 {
81 Size numberRates = taus_.size();
82
84 "B.rows()<> numberBumps_");
85
86 QL_REQUIRE(B.columns()==taus_.size(),
87 "B.columns()<> number of rates");
88
89 for (Size i =0; i < numberBumps_; ++i)
90 {
91 const Matrix& pseudo = pseudoBumped_[i];
92 driftsComputers_[i].compute(oldRates,
93 drifts_);
94
95 for (Size j =0; j < aliveIndex_; ++j)
96 B[i][j]=0.0;
97
98 for (Size j=aliveIndex_; j < numberRates; ++j)
99 {
100 bumpedRates_[j] = std::log(oldRates[j]+displacements_[j]);
101
102 for (Size k=0; k < factors_; ++k)
103 bumpedRates_[j] += -0.5*pseudo[j][k]*pseudo[j][k];
104
105 bumpedRates_[j] +=drifts_[j];
106
107 for (Size k=0; k < factors_; ++k)
108 bumpedRates_[j] += pseudo[j][k]*gaussians[k];
109
110 bumpedRates_[j] = std::exp(bumpedRates_[j]);
112 Real tmp = bumpedRates_[j] - newRates[j];
113
114 B[i][j] = tmp;
115 }
116 }
117
118 }
119
121 Size aliveIndex,
122 Size numeraire,
123 const std::vector<Time>& taus,
124 const std::vector<Matrix>& pseudoBumps,
125 std::vector<Spread> displacements)
126 : pseudoRoot_(pseudoRoot), aliveIndex_(aliveIndex), taus_(taus), pseudoBumps_(pseudoBumps),
127 displacements_(std::move(displacements)), numberBumps_(pseudoBumps.size()),
128 factors_(pseudoRoot.columns()),
129 // bumpedRates_(taus.size()),
130 e_(pseudoRoot.rows(), pseudoRoot.columns()), ratios_(taus_.size()) {
131 Size numberRates= taus.size();
132
133 QL_REQUIRE(aliveIndex == numeraire,
134 "we can do only do discretely compounding MM acount so aliveIndex must equal numeraire");
135
136 QL_REQUIRE(pseudoRoot_.rows()==numberRates,
137 "pseudoRoot_.rows()<> taus.size()");
138
139 QL_REQUIRE(displacements_.size()==numberRates,
140 "displacements_.size()<> taus.size()");
141
142
143 for (Size i=0; i < pseudoBumps.size(); ++i)
144 {
145 QL_REQUIRE(pseudoBumps[i].rows()==numberRates,
146 "pseudoBumps[i].rows()<> taus.size() with i =" << i);
147
148 QL_REQUIRE(pseudoBumps[i].columns()==factors_,
149 "pseudoBumps[i].columns()<> factors with i = " << i);
150
151
152
153
154 }
155
156 for (Size i=0; i < numberRates; ++i)
157 {
158 allDerivatives_.emplace_back(numberRates, factors_);
159 }
160 }
161
162
163 void RatePseudoRootJacobian::getBumps(const std::vector<Rate>& oldRates,
164 const std::vector<Real>& discountRatios,
165 const std::vector<Rate>& newRates,
166 const std::vector<Real>& gaussians,
167 Matrix& B)
168 {
169 Size numberRates = taus_.size();
170
171 QL_REQUIRE(B.rows() == numberBumps_, "we need B.rows() which is " << B.rows() << " to equal numberBumps_ which is " << numberBumps_);
172 QL_REQUIRE(B.columns() == numberRates, "we need B.columns() which is " << B.columns() << " to equal numberRates which is " << numberRates);
173
174
175 for (Size j=aliveIndex_; j < numberRates; ++j)
176 ratios_[j] = (oldRates[j] + displacements_[j])*discountRatios[j+1];
177
178 for (Size f=0; f < factors_; ++f)
179 {
180 e_[aliveIndex_][f] = 0;
181
182 for (Size j= aliveIndex_+1; j < numberRates; ++j)
183 e_[j][f] = e_[j-1][f] + ratios_[j-1]*pseudoRoot_[j-1][f];
184 }
185
186
187 for (Size f=0; f < factors_; ++f)
188 for (Size j=aliveIndex_; j < numberRates; ++j)
189 {
190 for (Size k= aliveIndex_; k < j ; ++k)
191 allDerivatives_[j][k][f] = newRates[j]*ratios_[k]*taus_[k]*pseudoRoot_[j][f];
192
193 // GG don't seem to have the 2, this term is miniscule in any case
194 Real tmp = //2*
195 2*ratios_[j]*taus_[j]*pseudoRoot_[j][f];
196 tmp -= pseudoRoot_[j][f];
197 tmp += e_[j][f]*taus_[j];
198 tmp += gaussians[f];
199 tmp *= (newRates[j]+displacements_[j]);
200
201
202 allDerivatives_[j][j][f] =tmp;
203
204 for (Size k= j+1; k < numberRates ; ++k)
205 allDerivatives_[j][k][f]=0.0;
206
207
208 }
209
210
211 for (Size i =0; i < numberBumps_; ++i)
212 {
213 Size j=0;
214
215 for (; j < aliveIndex_; ++j)
216 {
217 B[i][j]=0.0;
218 }
219 for (; j < numberRates; ++j)
220 {
221 Real sum =0.0;
222
223 for (Size k=aliveIndex_; k < numberRates; ++k)
224 for (Size f=0; f < factors_; ++f)
225 sum += pseudoBumps_[i][k][f]*allDerivatives_[j][k][f];
226 B[i][j] =sum;
227
228 }
229 }
230
231 }
232
234 const Matrix& pseudoRoot,
235 Size aliveIndex,
236 Size numeraire,
237 const std::vector<Time>& taus,
238 std::vector<Spread> displacements)
239 : pseudoRoot_(pseudoRoot), aliveIndex_(aliveIndex), taus_(taus),
240 displacements_(std::move(displacements)), factors_(pseudoRoot.columns()),
241 // bumpedRates_(taus.size()),
242 e_(pseudoRoot.rows(), pseudoRoot.columns()), ratios_(taus_.size()) {
243 Size numberRates= taus.size();
244
245 QL_REQUIRE(aliveIndex == numeraire,
246 "we can do only do discretely compounding MM acount so aliveIndex must equal numeraire");
247
248 QL_REQUIRE(pseudoRoot_.rows()==numberRates,
249 "pseudoRoot_.rows()<> taus.size()");
250
251 QL_REQUIRE(displacements_.size()==numberRates,
252 "displacements_.size()<> taus.size()");
253 }
254
255
256 void RatePseudoRootJacobianAllElements::getBumps(const std::vector<Rate>& oldRates,
257 const std::vector<Real>& discountRatios,
258 const std::vector<Rate>& newRates,
259 const std::vector<Real>& gaussians,
260 std::vector<Matrix>& B)
261 {
262 Size numberRates = taus_.size();
263
264 QL_REQUIRE(B.size() == numberRates, "we need B.size() which is " << B.size() << " to equal numberRates which is " << numberRates);
265 for (auto& j : B)
266 QL_REQUIRE(j.columns() == factors_ && j.rows() == numberRates,
267 "we need B[j].rows() which is "
268 << j.rows() << " to equal numberRates which is " << numberRates
269 << " and B[j].columns() which is " << j.columns()
270 << " to be equal to factors which is " << factors_);
271
272
273 for (Size j = aliveIndex_; j < numberRates; ++j)
274 ratios_[j] = (oldRates[j] + displacements_[j]) * discountRatios[j + 1];
275
276 for (Size f = 0; f < factors_; ++f) {
277 e_[aliveIndex_][f] = 0;
278
279 for (Size j = aliveIndex_ + 1; j < numberRates; ++j)
280 e_[j][f] = e_[j - 1][f] + ratios_[j - 1] * pseudoRoot_[j - 1][f];
281 }
282
283 // nullify B for rates that have already reset
284 for (Size j=0; j < aliveIndex_; ++j)
285 for (Size k=0; k < numberRates; ++k)
286 for (Size f=0; f < factors_; ++f)
287 B[j][k][f] =0.0;
288
289
290 for (Size f=0; f < factors_; ++f)
291 for (Size j=aliveIndex_; j < numberRates; ++j)
292 {
293 for (Size k= aliveIndex_; k < j ; ++k)
294 B[j][k][f] = newRates[j]*ratios_[k]*taus_[k]*pseudoRoot_[j][f];
295
296 Real tmp = 2*ratios_[j]*taus_[j]*pseudoRoot_[j][f];
297 tmp -= pseudoRoot_[j][f];
298 tmp += e_[j][f]*taus_[j];
299 tmp += gaussians[f];
300 tmp *= (newRates[j]+displacements_[j]);
301
302
303 B[j][j][f] =tmp;
304
305 for (Size k= 0; k < aliveIndex_ ; ++k)
306 B[j][k][f]=0.0;
307
308 for (Size k= j+1; k < numberRates ; ++k)
309 B[j][k][f]=0.0;
310
311
312 }
313 }
314
315
316}
317
318
319
320
Matrix used in linear algebra.
Definition: matrix.hpp:41
Size rows() const
Definition: matrix.hpp:504
Size columns() const
Definition: matrix.hpp:508
Matrix pseudoRoot_
this data does not change after construction
RatePseudoRootJacobianAllElements(const Matrix &pseudoRoot, Size aliveIndex, Size numeraire, const std::vector< Time > &taus, std::vector< Spread > displacements)
void getBumps(const std::vector< Rate > &oldRates, const std::vector< Real > &oneStepDFs, const std::vector< Rate > &newRates, const std::vector< Real > &gaussians, std::vector< Matrix > &B)
Matrix pseudoRoot_
this data does not change after construction
std::vector< Matrix > allDerivatives_
workspace variables
void getBumps(const std::vector< Rate > &oldRates, const std::vector< Real > &oneStepDFs, const std::vector< Rate > &newRates, const std::vector< Real > &gaussians, Matrix &B)
RatePseudoRootJacobian(const Matrix &pseudoRoot, Size aliveIndex, Size numeraire, const std::vector< Time > &taus, const std::vector< Matrix > &pseudoBumps, std::vector< Spread > displacements)
Matrix pseudoRoot_
this data is always the same
RatePseudoRootJacobianNumerical(const Matrix &pseudoRoot, Size aliveIndex, Size numeraire, const std::vector< Time > &taus, const std::vector< Matrix > &pseudoBumps, const std::vector< Spread > &displacements)
void getBumps(const std::vector< Rate > &oldRates, const std::vector< Real > &oneStepDFs, const std::vector< Rate > &newRates, const std::vector< Real > &gaussians, Matrix &B)
std::vector< Real > drifts_
workspace variables
std::vector< LMMDriftCalculator > driftsComputers_
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
Definition: errors.hpp:117
QL_REAL Real
real number
Definition: types.hpp:50
std::size_t Size
size of a container
Definition: types.hpp:58
Definition: any.hpp:35
STL namespace.