QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
alphafinder.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) 2007 Mark Joshi
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/math/quadratic.hpp>
22#include <utility>
23
24namespace QuantLib {
25
26namespace
27{
28
29 template<class T>
30 Real Bisection(Real target,
31 Real low,
32 Real high,
33 Real tolerance,
34 T& theObject,
35 Real (T::*Value)(Real)) {
36
37 Real x=0.5*(low+high);
38 Real y=(theObject.*Value)(x);
39
40 do {
41 if (y < target) low = x;
42 else if (y > target) high = x;
43
44 x = 0.5*(low+high);
45 y = (theObject.*Value)(x);
46 } while ((std::fabs(high-low) > tolerance));
47
48 return x;
49 }
50
51 template<class T>
52 Real FindHighestOK(Real low,
53 Real high,
54 Real tolerance,
55 T& theObject,
56 bool (T::*Value)(Real)) {
57
58 Real x=0.5*(low+high);
59 bool ok=(theObject.*Value)(x);
60
61 do {
62 if (ok) low = x;
63 else high = x;
64
65 x = 0.5*(low+high);
66 ok = (theObject.*Value)(x);
67 } while ((std::fabs(high-low) > tolerance));
68
69 return x;
70 }
71
72 template<class T>
73 Real FindLowestOK(Real low,
74 Real high,
75 Real tolerance,
76 T& theObject,
77 bool (T::*Value)(Real)) {
78
79 Real x=0.5*(low+high);
80 bool ok=(theObject.*Value)(x);
81
82 do {
83 if (ok) high = x;
84 else low = x;
85
86 x = 0.5*(low+high);
87 ok = (theObject.*Value)(x);
88 } while ( (std::fabs(high-low) > tolerance) );
89
90 return x;
91 }
92
93
94 template<class T>
95 Real Minimize(Real low,
96 Real high,
97 Real tolerance,
98 T& theObject,
99 Real (T::*Value)(Real),
100 bool (T::*Condition)(Real),
101 bool& failed) {
102
103 Real leftValue = (theObject.*Value)(low);
104 Real rightValue = (theObject.*Value)(high);
105 Real W = 0.5*(3.0-std::sqrt(5.0));
106 Real x=W*low+(1-W)*high;
107 Real midValue = (theObject.*Value)(x);
108
109 failed = true;
110
111 while(high - low > tolerance) {
112
113 if (x - low > high -x) // left interval is bigger
114 {
115 Real tentativeNewMid = W*low+(1-W)*x;
116 Real tentativeNewMidValue = (theObject.*Value)(tentativeNewMid);
117 bool conditioner = (theObject.*Condition)(tentativeNewMidValue);
118 if (!conditioner) {
119 if ((theObject.*Condition)(x))
120 return x;
121 else
122 if (leftValue < rightValue)
123 return low;
124 else
125 return high;
126 }
127
128 if (tentativeNewMidValue < midValue) // go left
129 {
130 high =x;
131 rightValue = midValue;
132 x = tentativeNewMid;
133 midValue = tentativeNewMidValue;
134 }
135 else // go right
136 {
137 low = tentativeNewMid;
138 leftValue = tentativeNewMidValue;
139 }
140 }
141 else
142 {
143 Real tentativeNewMid = W*x+(1-W)*high;
144 Real tentativeNewMidValue = (theObject.*Value)(tentativeNewMid);
145 bool conditioner = (theObject.*Condition)(tentativeNewMidValue);
146 if (!conditioner) {
147 if ((theObject.*Condition)(x))
148 return x;
149 else
150 if (leftValue < rightValue)
151 return low;
152 else
153 return high;
154 }
155
156 if (tentativeNewMidValue < midValue) // go right
157 {
158 low =x;
159 leftValue = midValue;
160 x = tentativeNewMid;
161 midValue = tentativeNewMidValue;
162 }
163 else // go left
164 {
165 high = tentativeNewMid;
166 rightValue = tentativeNewMidValue;
167 }
168 }
169
170
171
172
173 }
174 failed = false;
175 return x;
176 }
177}
178
179AlphaFinder::AlphaFinder(ext::shared_ptr<AlphaForm> parametricform)
180: parametricform_(std::move(parametricform)) {}
181
182
184 Real cov = 0.0;
185 parametricform_->setAlpha(alpha);
186
187 for (Integer i = 0; i < stepindex_ + 1; ++i) {
188 Real vol1 = ratetwohomogeneousvols_[i] * (*parametricform_)(i);
189 cov += vol1 * rateonevols_[i] * correlations_[i];
190 }
191
192 cov *= 2 * w0_ * w1_;
193 return cov;
194 }
195
196
198 Real var =0.0;
199 parametricform_->setAlpha(alpha);
200
201 for (Integer i=0; i < stepindex_+1; ++i) {
202 Real vol = ratetwohomogeneousvols_[i]*(*parametricform_)(i);
203 var+= vol*vol;
204 }
205
206 var *= w1_*w1_;
207 return var;
208 }
209
211 Real dum1, dum2, dum3;
218 dum1,
219 dum2,
220 dum3,
222
223 Real result=0.0;
224 for (Size i=0; i<=static_cast<Size>(stepindex_)+1; ++i) {
226 result +=val*val;
227 }
228
229 return result;
230}
231
233 Integer stepindex,
234 const std::vector<Volatility>& ratetwohomogeneousvols,
235 Real quadraticPart,
236 Real linearPart,
237 Real constantPart,
238 Real& alpha,
239 Real& a,
240 Real& b,
241 std::vector<Volatility>& ratetwovols) {
242 alpha = alphaFound;
243 quadratic q2(quadraticPart, linearPart, constantPart-targetVariance_ );
244 parametricform_->setAlpha(alpha);
245 Real y; // dummy
246 q2.roots(a,y);
247
248 Real varSoFar=0.0;
249 for (Integer i =0; i < stepindex+1; ++i) {
250 ratetwovols[i] = ratetwohomogeneousvols[i] *
251 (*parametricform_)(i) * a;
252 varSoFar += ratetwovols[i]* ratetwovols[i];
253 }
254
255 Real VarToFind = totalVar_-varSoFar;
256 if (VarToFind < 0)
257 return false;
258 Real requiredSd = std::sqrt(VarToFind);
259 b = requiredSd / (ratetwohomogeneousvols[stepindex+1] *
260 (*parametricform_)(stepindex));
261 ratetwovols[stepindex+1] = requiredSd;
262 return true;
263 }
264
266
269
271 Real valueAtTP =q.valueAtTurningPoint();
272 return valueAtTP;
273 }
274
276 return -valueAtTurningPoint(alpha);
277 }
278
281 if (!aExists)
282 return false;
283
284 Real dum1, dum2, dum3;
285 return finalPart(alpha,
291 dum1,
292 dum2,
293 dum3,
295 }
296
298 Integer stepindex, // caplet index
299 const std::vector<Volatility>& rateonevols,
300 const std::vector<Volatility>& ratetwohomogeneousvols,
301 const std::vector<Real>& correlations,
302 Real w0,
303 Real w1,
304 Real targetVariance,
305 Real tolerance,
306 Real alphaMax,
307 Real alphaMin,
308 Integer steps,
309 Real& alpha,
310 Real& a,
311 Real& b,
312 std::vector<Volatility>& ratetwovols) {
313 stepindex_=stepindex;
314 rateonevols_=rateonevols;
315 ratetwohomogeneousvols_=ratetwohomogeneousvols;
316 correlations_=correlations;
317 w0_=w0;
318 w1_=w1;
319 totalVar_=0;
320 for (Size i=0; i <=static_cast<Size>(stepindex)+1; ++i)
321 totalVar_+=ratetwohomogeneousvols[i]*ratetwohomogeneousvols[i];
322 targetVariance_ = targetVariance;
323
324 // constant part will not depend on alpha
325
326 constantPart_ =0.0;
327 for (Integer i=0; i < stepindex+1; ++i)
328 constantPart_+=rateonevols[i]*rateonevols[i];
329
330 constantPart_ *= w0*w0;
331
332 // compute linear part with initial alpha
333 Real valueAtTP = valueAtTurningPoint(alpha0);
334
335 if (valueAtTP <= targetVariance) {
336 finalPart(alpha0,
337 stepindex,
338 ratetwohomogeneousvols,
342 alpha,
343 a,
344 b,
345 ratetwovols);
346 return true;
347 }
348
349 // we now have to solve
350 Real bottomValue = valueAtTurningPoint(alphaMin);
351 Real bottomAlpha = alphaMin;
352 Real topValue = valueAtTurningPoint(alphaMax);
353 Real topAlpha = alphaMax;
354 Real bilimit = alpha0;
355
356 if (bottomValue > targetVariance && topValue > targetVariance) {
357 // see if if ok at some intermediate point by stepping through
358 Integer i=1;
359 while ( i < steps && topValue> targetVariance) {
360 topAlpha = alpha0 + (alphaMax-alpha0)*(i+0.0)/(steps+0.0);
361 topValue=valueAtTurningPoint(topAlpha);
362 ++i;
363 }
364
365 if (topValue <= targetVariance)
366 bilimit = alpha0 + (topAlpha-alpha0)*(i-2.0)/(steps+0.0);
367 }
368
369 if (bottomValue > targetVariance && topValue > targetVariance) {
370 // see if if ok at some intermediate point by stepping through
371 Integer i=1;
372 while ( i < steps && topValue> targetVariance) {
373 bottomAlpha = alpha0 + (alphaMin-alpha0)*(i+0.0)/(steps+0.0);
374 bottomValue=valueAtTurningPoint(bottomAlpha);
375 ++i;
376 }
377
378 if (bottomValue <= targetVariance)
379 bilimit = alpha0 +(bottomAlpha-alpha0)*(i-2.0)/(steps+0.0);
380 }
381
382 if (bottomValue > targetVariance && topValue > targetVariance)
383 return false;
384
385 if (bottomValue <= targetVariance) {
386 // then find root of increasing function
387 // (or as if increasing function)
389 targetVariance,
390 bottomAlpha,
391 bilimit,
392 tolerance,
393 *this,
395 } else {
396 // find root of decreasing function (or as if decreasing function)
398 -targetVariance,
399 bilimit,
400 topAlpha,
401 tolerance,
402 *this,
404 }
406 stepindex,
407 ratetwohomogeneousvols,
411 alpha,
412 a,
413 b,
414 ratetwovols);
415 return true;
416 }
417
419 Real alpha0,
420 Integer stepindex, // caplet index
421 const std::vector<Volatility>& rateonevols,
422 const std::vector<Volatility>& ratetwohomogeneousvols,
423 const std::vector<Real>& correlations,
424 Real w0,
425 Real w1,
426 Real targetVariance,
427 Real tolerance,
428 Real alphaMax,
429 Real alphaMin,
430 Integer steps,
431 Real& alpha,
432 Real& a,
433 Real& b,
434 std::vector<Volatility>& ratetwovols) {
435
436 stepindex_=stepindex;
437 rateonevols_=rateonevols;
438 ratetwohomogeneousvols_=ratetwohomogeneousvols;
440 correlations_=correlations;
441 w0_=w0;
442 w1_=w1;
443 totalVar_=0;
444 for (Size i=0; i <=static_cast<Size>(stepindex)+1; ++i)
445 totalVar_+=ratetwohomogeneousvols[i]*ratetwohomogeneousvols[i];
446 targetVariance_=targetVariance;
447
448 // constant part will not depend on alpha
449
450 constantPart_ =0.0;
451 for (Integer i=0; i < stepindex+1; ++i)
452 constantPart_+=rateonevols[i]*rateonevols[i];
453
454 constantPart_ *= w0*w0;
455
456 Real alpha1 = alphaMin;
457 Real alpha2 = alphaMax;
458
459 // compute linear part with initial alpha
460 bool alpha0OK = testIfSolutionExists(alpha0);
461 bool alphaMaxOK = testIfSolutionExists(alphaMax);
462 bool alphaMinOK = testIfSolutionExists(alphaMin);
463
464 bool foundOKPoint = alpha0OK || alphaMaxOK || alphaMinOK;
465
466 if (foundOKPoint) {
467 if (!alphaMinOK) {
468 // lower alpha is bad
469 if (alpha0OK) {
470 // must die somewhere in between
471 alpha1 = FindLowestOK<AlphaFinder>(
472 alphaMin,
473 alpha0,
474 tolerance,
475 *this,
477 } else {
478 // alphaMaxOK must be true to get here
479 alpha1 = FindLowestOK<AlphaFinder>(
480 alpha0,
481 alphaMax,
482 tolerance,
483 *this,
485 }
486 }
487
488
489 if (!alphaMaxOK) {
490 // higher alpha is bad
491 alpha2 = FindHighestOK<AlphaFinder>(
492 alpha1,
493 alphaMax,
494 tolerance,
495 *this,
497 } else
498 alpha2= alphaMax;
499 }
500 else {
501 // ok let's see if we can find a value of alpha that works
502 bool foundUpOK = false;
503 bool foundDownOK = false;
504 Real alphaUp = alpha0, alphaDown = alpha0;
505 Real stepSize = (alphaMax-alpha0)/steps;
506
507 for (Size j=0;
508 j<static_cast<Size>(steps) && !foundUpOK && !foundDownOK;
509 ++j) {
510 alphaUp = alpha0+j*stepSize;
511 foundUpOK=testIfSolutionExists(alphaUp);
512 alphaDown = alpha0-j*stepSize;
513 foundDownOK=testIfSolutionExists(alphaDown);
514 }
515 foundOKPoint = foundUpOK || foundDownOK;
516 if (!foundOKPoint)
517 return false;
518
519 if (foundUpOK) {
520 alpha1 = alphaUp;
521 alpha2 = FindHighestOK<AlphaFinder>(
522 alpha1,
523 alphaMax,
524 tolerance,
525 *this,
527 } else {
528 alpha2 = alphaDown;
529 alpha1 = FindLowestOK<AlphaFinder>(
530 alphaMin,
531 alpha2,
532 tolerance,
533 *this,
535 }
536 }
537
538 // we have now found alpha1, alpha2 such that solution exists
539 // at endpoints. we now want to minimize within that interval
540 bool failed;
541 alpha = Minimize<AlphaFinder>(
542 alpha1,
543 alpha2,
544 tolerance,
545 *this,
548 failed) ;
549
551 stepindex,
552 ratetwohomogeneousvols,
556 alpha,
557 a,
558 b,
559 ratetwovols);
560
561 return true;;
562 }
563}
Real minusValueAtTurningPoint(Real alpha)
Real homogeneityfailure(Real alpha)
std::vector< Volatility > ratetwohomogeneousvols_
Definition: alphafinder.hpp:86
bool testIfSolutionExists(Real alpha)
Real valueAtTurningPoint(Real alpha)
bool solveWithMaxHomogeneity(Real alpha0, Integer stepindex, const std::vector< Volatility > &rateonevols, const std::vector< Volatility > &ratetwohomogeneousvols, const std::vector< Real > &correlations, Real w0, Real w1, Real targetVariance, Real tolerance, Real alphaMax, Real alphaMin, Integer steps, Real &alpha, Real &a, Real &b, std::vector< Volatility > &ratetwovols)
bool finalPart(Real alphaFound, Integer stepindex, const std::vector< Volatility > &ratetwohomogeneousvols, Real quadraticPart, Real linearPart, Real constantPart, Real &alpha, Real &a, Real &b, std::vector< Volatility > &ratetwovols)
bool solve(Real alpha0, Integer stepindex, const std::vector< Volatility > &rateonevols, const std::vector< Volatility > &ratetwohomogeneousvols, const std::vector< Real > &correlations, Real w0, Real w1, Real targetVariance, Real tolerance, Real alphaMax, Real alphaMin, Integer steps, Real &alpha, Real &a, Real &b, std::vector< Volatility > &ratetwovols)
ext::shared_ptr< AlphaForm > parametricform_
Definition: alphafinder.hpp:84
Real computeQuadraticPart(Real alpha)
std::vector< Volatility > rateonevols_
Definition: alphafinder.hpp:86
std::vector< Real > correlations_
Definition: alphafinder.hpp:88
Real computeLinearPart(Real alpha)
std::vector< Volatility > putativevols_
Definition: alphafinder.hpp:87
AlphaFinder(ext::shared_ptr< AlphaForm > parametricform)
Bisection 1-D solver
Definition: bisection.hpp:37
bool roots(Real &x, Real &y) const
Definition: quadratic.cpp:44
ext::function< Real(Real)> b
QL_REAL Real
real number
Definition: types.hpp:50
QL_INTEGER Integer
integer number
Definition: types.hpp:35
std::size_t Size
size of a container
Definition: types.hpp:58
Definition: any.hpp:35
STL namespace.
ext::shared_ptr< YieldTermStructure > q
quadratic formula
Real alpha
Definition: sabr.cpp:200