QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
pathwiseaccountingengine.cpp
Go to the documentation of this file.
1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4Copyright (C) 2008 Mark Joshi
5
6This file is part of QuantLib, a free-software/open-source library
7for financial quantitative analysts and developers - http://quantlib.org/
8
9QuantLib is free software: you can redistribute it and/or modify it
10under the terms of the QuantLib license. You should have received a
11copy 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
15This program is distributed in the hope that it will be useful, but WITHOUT
16ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17FOR A PARTICULAR PURPOSE. See the license for more details.
18*/
19
26#include <algorithm>
27#include <utility>
28
29namespace QuantLib {
30
32 ext::shared_ptr<LogNormalFwdRateEuler> evolver, // method relies heavily on LMM Euler
34 ext::shared_ptr<MarketModel> pseudoRootStructure, // we need pseudo-roots and displacements
35 Real initialNumeraireValue)
36 : evolver_(std::move(evolver)), product_(product),
37 pseudoRootStructure_(std::move(pseudoRootStructure)),
38 initialNumeraireValue_(initialNumeraireValue), numberProducts_(product->numberOfProducts()),
39 doDeflation_(!product->alreadyDeflated()), numerairesHeld_(product->numberOfProducts()),
40 numberCashFlowsThisStep_(product->numberOfProducts()),
41 cashFlowsGenerated_(product->numberOfProducts()),
42 deflatorAndDerivatives_(pseudoRootStructure_->numberOfRates() + 1) {
43
44 numberRates_ = pseudoRootStructure_->numberOfRates();
45 numberSteps_ = pseudoRootStructure_->numberOfSteps();
46
48
49
50
52
53 for (Size i=0; i <= numberSteps_; ++i)
54 Discounts_[i][0] = 1.0;
55
56
57 V_.reserve(numberProducts_);
58
59 Matrix modelCashFlowIndex(product_->possibleCashFlowTimes().size(), numberRates_+1);
60
61
63
64 for (Size i=0; i<numberProducts_; ++i)
65 {
66 cashFlowsGenerated_[i].resize(
67 product_->maxNumberOfCashFlowsPerProductPerStep());
68
69 for (auto& j : cashFlowsGenerated_[i])
70 j.amount.resize(numberRates_ + 1);
71
72 numberCashFlowsThisIndex_[i].resize(product_->possibleCashFlowTimes().size());
73
74 V_.push_back(VModel);
75
76
77 totalCashFlowsThisIndex_.push_back(modelCashFlowIndex);
78 }
79
80 LIBORRatios_ = VModel;
82 LIBORRates_ =VModel;
83
84
85
86
87 const std::vector<Time>& cashFlowTimes =
88 product_->possibleCashFlowTimes();
89 numberCashFlowTimes_ = cashFlowTimes.size();
90
91 const std::vector<Time>& rateTimes = product_->evolution().rateTimes();
92 const std::vector<Time>& evolutionTimes = product_->evolution().evolutionTimes();
93 discounters_.reserve(cashFlowTimes.size());
94
95 for (Real cashFlowTime : cashFlowTimes)
96 discounters_.emplace_back(cashFlowTime, rateTimes);
97
98
99 // need to check that we are in money market measure
100
101
102 // we need to allocate cash-flow times to steps, i.e. what is the last step completed before a flow occurs
103 // what we really need is for each step, what cash flow time indices to look at
104
106
107 for (Size i=0; i < numberCashFlowTimes_; ++i)
108 {
109 auto it =
110 std::upper_bound(evolutionTimes.begin(), evolutionTimes.end(), cashFlowTimes[i]);
111 if (it != evolutionTimes.begin())
112 --it;
113 Size index = it - evolutionTimes.begin();
114 cashFlowIndicesThisStep_[index].push_back(i);
115 }
116
118 }
119
121 {
122
123 const std::vector<Real> initialForwards_(pseudoRootStructure_->initialRates());
124 currentForwards_ = initialForwards_;
125 // clear accumulation variables
126 for (Size i=0; i < numberProducts_; ++i)
127 {
128 numerairesHeld_[i]=0.0;
129
130 for (Size j=0; j < numberCashFlowTimes_; ++j)
131 {
133
134 for (Size k=0; k <= numberRates_; ++k)
135 totalCashFlowsThisIndex_[i][j][k] =0.0;
136 }
137
138 for (Size l=0; l< numberRates_; ++l)
139 for (Size m=0; m <= numberSteps_; ++m)
140 V_[i][m][l] =0.0;
141
142 }
143
144
145
146 Real weight = evolver_->startNewPath();
147 product_->reset();
148
149 Size thisStep;
150
151 bool done = false;
152 do {
153 thisStep = evolver_->currentStep();
154 Size storeStep = thisStep+1;
155 weight *= evolver_->advanceStep();
156
157 done = product_->nextTimeStep(evolver_->currentState(),
160
162 currentForwards_ = evolver_->currentState().forwardRates();
163
164 for (unsigned long i=0; i < numberRates_; ++i)
165 {
166 Real x= evolver_->currentState().discountRatio(i+1,i);
167 StepsDiscountsSquared_[storeStep][i] = x*x;
168
169 LIBORRatios_[storeStep][i] = currentForwards_[i]/lastForwards_[i];
170 LIBORRates_[storeStep][i] = currentForwards_[i];
171 Discounts_[storeStep][i+1] = evolver_->currentState().discountRatio(i+1,0);
172 }
173
174 // for each product...
175 for (Size i=0; i<numberProducts_; ++i)
176 {
177 // ...and each cash flow...
178 for (Size j=0; j<numberCashFlowsThisStep_[i]; ++j)
179 {
180 Size k = cashFlowsGenerated_[i][j].timeIndex;
182
183 for (Size l=0; l <= numberRates_; ++l)
184 totalCashFlowsThisIndex_[i][k][l] += cashFlowsGenerated_[i][j].amount[l]*weight;
185
186 }
187 }
188
189
190 } while (!done);
191
192 // ok we've gathered cash-flows, still have to backwards computation
193
194 Size factors = pseudoRootStructure_->numberOfFactors();
195 const std::vector<Time>& taus= pseudoRootStructure_->evolution(). rateTaus();
196
197 bool flowsFound = false;
198
199 Integer finalStepDone = thisStep;
200
201 for (Integer currentStep = numberSteps_-1; currentStep >=0 ; --currentStep) // must be a signed type as we go negative
202 {
203 Integer stepToUse = std::min<Integer>(currentStep, finalStepDone)+1;
204
205 for (Size k=0; k < cashFlowIndicesThisStep_[currentStep].size(); ++k)
206 {
207 Size cashFlowIndex =cashFlowIndicesThisStep_[currentStep][k];
208
209 // first check to see if anything actually happened before spending time on computing stuff
210 bool noFlows = true;
211 for (Size l=0; l < numberProducts_ && noFlows; ++l)
212 noFlows = noFlows && (numberCashFlowsThisIndex_[l][cashFlowIndex] ==0);
213
214 flowsFound = flowsFound || !noFlows;
215
216 if (!noFlows)
217 {
218 if (doDeflation_)
219 discounters_[cashFlowIndex].getFactors(LIBORRates_, Discounts_,stepToUse, deflatorAndDerivatives_); // get amount to discount cash flow by and amount to multiply its derivatives by
220
221 for (Size j=0; j < numberProducts_; ++j)
222 {
223 if (numberCashFlowsThisIndex_[j][cashFlowIndex] > 0)
224 {
225 Real deflatedCashFlow = totalCashFlowsThisIndex_[j][cashFlowIndex][0];
226 if (doDeflation_)
227 deflatedCashFlow *= deflatorAndDerivatives_[0];
228 //cashFlowsGenerated_[j][cashFlowIndex].amount[0]*deflatorAndDerivatives_[0];
229 numerairesHeld_[j] += deflatedCashFlow;
230
231 for (Size i=1; i <= numberRates_; ++i)
232 {
233 Real thisDerivative = totalCashFlowsThisIndex_[j][cashFlowIndex][i];
234 if (doDeflation_)
235 {
236 thisDerivative *= deflatorAndDerivatives_[0];
237 thisDerivative += totalCashFlowsThisIndex_[j][cashFlowIndex][0]*deflatorAndDerivatives_[i];
238 }
239
240 V_[j][stepToUse][i-1] += thisDerivative; // zeroth row of V is t =0 not t_0
241 }
242 }
243 }
244 }
245 }
246
247 // need to do backwards updating
248 if (flowsFound)
249 {
250 Integer nextStepToUse = std::min<Integer>(currentStep-1, finalStepDone);
251 Integer nextStepIndex = nextStepToUse+1;
252 if (nextStepIndex != stepToUse) // then we need to update V
253 {
254
255 const Matrix& thisPseudoRoot_= pseudoRootStructure_->pseudoRoot(currentStep);
256
257 for (Size i=0; i < numberProducts_; ++i)
258 {
259 // compute partials
260 for (Size f=0; f < factors; ++f)
261 {
262 Real libor = LIBORRates_[stepToUse][numberRates_-1];
263 Real V = V_[i][stepToUse][numberRates_-1];
264 Real pseudo = thisPseudoRoot_[numberRates_-1][f];
265 Real thisPartialTerm = libor*V*pseudo;
266 partials_[f][numberRates_-1] = thisPartialTerm;
267
268 for (Integer r = numberRates_-2; r >=0 ; --r)
269 {
270 Real thisPartialTermr = LIBORRates_[stepToUse][r]*V_[i][stepToUse][r]*thisPseudoRoot_[r][f];
271
272 partials_[f][r] = partials_[f][r+1] + thisPartialTermr;
273
274 }
275 }
276 for (Size j=0; j < numberRates_; ++j)
277 {
278 Real nextV = V_[i][stepToUse][j] * LIBORRatios_[stepToUse][j];
279 V_[i][nextStepIndex][j] = nextV;
280
281 Real summandTerm = 0.0;
282 for (Size f=0; f < factors; ++f)
283 summandTerm += thisPseudoRoot_[j][f]*partials_[f][j];
284
285 summandTerm *= taus[j]*StepsDiscountsSquared_[stepToUse][j];
286
287 V_[i][nextStepIndex][j] += summandTerm;
288
289 }
290 }
291
292 }
293 }
294
295
296
297
298 }
299
300 // write answer into values
301
302 for (Size i=0; i < numberProducts_; ++i)
303 {
305 for (Size j=0; j < numberRates_; ++j)
306 values[(i+1)*numberProducts_+j] = V_[i][0][j]*initialNumeraireValue_;
307 }
308
309 return 1.0; // we have put the weight in already, this results in lower variance since weight changes along the path
310 }
311
313 Size numberOfPaths)
314 {
315 std::vector<Real> values(product_->numberOfProducts()*(numberRates_+1));
316 for (Size i=0; i<numberOfPaths; ++i)
317 {
318 Real weight = singlePathValues(values);
319 stats.add(values,weight);
320 }
321 }
322
323////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
324
325////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
326
328 ext::shared_ptr<LogNormalFwdRateEuler> evolver, // method relies heavily on LMM Euler
330 ext::shared_ptr<MarketModel> pseudoRootStructure, // we need pseudo-roots and displacements
331 const std::vector<std::vector<Matrix> >& vegaBumps,
332 Real initialNumeraireValue)
333 : evolver_(std::move(evolver)), product_(product),
334 pseudoRootStructure_(std::move(pseudoRootStructure)),
335 initialNumeraireValue_(initialNumeraireValue), numberProducts_(product->numberOfProducts()),
336 doDeflation_(!product->alreadyDeflated()), numerairesHeld_(product->numberOfProducts()),
337 numberCashFlowsThisStep_(product->numberOfProducts()),
338 cashFlowsGenerated_(product->numberOfProducts()),
339 stepsDiscounts_(pseudoRootStructure_->numberOfRates() + 1),
340 vegasThisPath_(product->numberOfProducts(), vegaBumps[0].size()),
341 deflatorAndDerivatives_(pseudoRootStructure_->numberOfRates() + 1) {
342
343 stepsDiscounts_[0]=1.0;
344
345 numberRates_ = pseudoRootStructure_->numberOfRates();
346 numberSteps_ = pseudoRootStructure_->numberOfSteps();
348
349
350 const EvolutionDescription& evolution = pseudoRootStructure_->evolution();
351 numeraires_ = moneyMarketMeasure(evolution);
352
353
354 QL_REQUIRE(vegaBumps.size() == numberSteps_, "we need one vector of vega bumps for each step.");
355
356 numberBumps_ = vegaBumps[0].size();
357
358 for (Size i =0; i < numberSteps_; ++i)
359 {
360 Size thisSize = vegaBumps[i].size();
361 QL_REQUIRE(thisSize == numberBumps_,"We must have precisely the same number of bumps for each step.");
362 jacobianComputers_.emplace_back(
363 pseudoRootStructure_->pseudoRoot(i), evolution.firstAliveRate()[i], numeraires_[i],
364 evolution.rateTaus(), vegaBumps[i], pseudoRootStructure_->displacements());
365
366 jacobiansThisPaths_.emplace_back(numberBumps_, pseudoRootStructure_->numberOfRates());
367 }
368
369
370
372
373
374
376
377 for (Size i=0; i <= numberSteps_; ++i)
378 Discounts_[i][0] = 1.0;
379
380
381 V_.reserve(numberProducts_);
382
383 Matrix modelCashFlowIndex(product_->possibleCashFlowTimes().size(), numberRates_+1);
384
385
387
388 for (Size i=0; i<numberProducts_; ++i)
389 {
390 cashFlowsGenerated_[i].resize(
391 product_->maxNumberOfCashFlowsPerProductPerStep());
392
393 for (auto& j : cashFlowsGenerated_[i])
394 j.amount.resize(numberRates_ + 1);
395
396 numberCashFlowsThisIndex_[i].resize(product_->possibleCashFlowTimes().size());
397
398 V_.push_back(VModel);
399
400
401 totalCashFlowsThisIndex_.push_back(modelCashFlowIndex);
402 }
403
404 LIBORRatios_ = VModel;
405 StepsDiscountsSquared_ = VModel;
406 LIBORRates_ =VModel;
407
408
409
410
411 const std::vector<Time>& cashFlowTimes =
412 product_->possibleCashFlowTimes();
413 numberCashFlowTimes_ = cashFlowTimes.size();
414
415 const std::vector<Time>& rateTimes = product_->evolution().rateTimes();
416 const std::vector<Time>& evolutionTimes = product_->evolution().evolutionTimes();
417 discounters_.reserve(cashFlowTimes.size());
418
419 for (Real cashFlowTime : cashFlowTimes)
420 discounters_.emplace_back(cashFlowTime, rateTimes);
421
422
423 // need to check that we are in money market measure
424
425
426 // we need to allocate cash-flow times to steps, i.e. what is the last step completed before a flow occurs
427 // what we really need is for each step, what cash flow time indices to look at
428
430
431 for (Size i=0; i < numberCashFlowTimes_; ++i)
432 {
433 auto it =
434 std::upper_bound(evolutionTimes.begin(), evolutionTimes.end(), cashFlowTimes[i]);
435 if (it != evolutionTimes.begin())
436 --it;
437 Size index = it - evolutionTimes.begin();
438 cashFlowIndicesThisStep_[index].push_back(i);
439 }
440
442 }
443
445 {
446
447 const std::vector<Real>& initialForwards_(pseudoRootStructure_->initialRates());
448 currentForwards_ = initialForwards_;
449 // clear accumulation variables
450 for (Size i=0; i < numberProducts_; ++i)
451 {
452 numerairesHeld_[i]=0.0;
453
454 for (Size j=0; j < numberCashFlowTimes_; ++j)
455 {
457
458 for (Size k=0; k <= numberRates_; ++k)
459 totalCashFlowsThisIndex_[i][j][k] =0.0;
460 }
461
462 for (Size l=0; l< numberRates_; ++l)
463 for (Size m=0; m <= numberSteps_; ++m)
464 V_[i][m][l] =0.0;
465
466 for (Size p=0; p < numberBumps_; ++p)
467 vegasThisPath_[i][p] =0.0;
468
469 }
470
471
472
473 Real weight = evolver_->startNewPath();
474 product_->reset();
475
476 Size thisStep;
477
478 bool done = false;
479 do {
480 thisStep = evolver_->currentStep();
481 Size storeStep = thisStep+1;
482 weight *= evolver_->advanceStep();
483
484 done = product_->nextTimeStep(evolver_->currentState(),
487
489 currentForwards_ = evolver_->currentState().forwardRates();
490
491 for (unsigned long i=0; i < numberRates_; ++i)
492 {
493 Real x= evolver_->currentState().discountRatio(i+1,i);
494 stepsDiscounts_[i+1] = x;
495 StepsDiscountsSquared_[storeStep][i] = x*x;
496
497 LIBORRatios_[storeStep][i] = currentForwards_[i]/lastForwards_[i];
498 LIBORRates_[storeStep][i] = currentForwards_[i];
499 Discounts_[storeStep][i+1] = evolver_->currentState().discountRatio(i+1,0);
500 }
501
502 jacobianComputers_[thisStep].getBumps(lastForwards_,
505 evolver_->browniansThisStep(),
506 jacobiansThisPaths_[thisStep]);
507
508
509
510 // for each product...
511 for (Size i=0; i<numberProducts_; ++i)
512 {
513 // ...and each cash flow...
514 for (Size j=0; j<numberCashFlowsThisStep_[i]; ++j)
515 {
516 Size k = cashFlowsGenerated_[i][j].timeIndex;
518
519 for (Size l=0; l <= numberRates_; ++l)
520 totalCashFlowsThisIndex_[i][k][l] += cashFlowsGenerated_[i][j].amount[l]*weight;
521
522 }
523 }
524
525
526 } while (!done);
527
528 // ok we've gathered cash-flows, still have to backwards computation
529
530 Size factors = pseudoRootStructure_->numberOfFactors();
531 const std::vector<Time>& taus= pseudoRootStructure_->evolution(). rateTaus();
532
533 bool flowsFound = false;
534
535 Integer finalStepDone = thisStep;
536
537 for (Integer currentStep = numberSteps_-1; currentStep >=0 ; --currentStep) // must be a signed type as we go negative
538 {
539 Integer stepToUse = std::min<Integer>(currentStep, finalStepDone)+1;
540
541 for (Size k=0; k < cashFlowIndicesThisStep_[currentStep].size(); ++k)
542 {
543 Size cashFlowIndex =cashFlowIndicesThisStep_[currentStep][k];
544
545 // first check to see if anything actually happened before spending time on computing stuff
546 bool noFlows = true;
547 for (Size l=0; l < numberProducts_ && noFlows; ++l)
548 noFlows = noFlows && (numberCashFlowsThisIndex_[l][cashFlowIndex] ==0);
549
550 flowsFound = flowsFound || !noFlows;
551
552 if (!noFlows)
553 {
554 if (doDeflation_)
555 discounters_[cashFlowIndex].getFactors(LIBORRates_, Discounts_,stepToUse, deflatorAndDerivatives_); // get amount to discount cash flow by and amount to multiply its derivatives by
556
557 for (Size j=0; j < numberProducts_; ++j)
558 {
559 if (numberCashFlowsThisIndex_[j][cashFlowIndex] > 0)
560 {
561 Real deflatedCashFlow = totalCashFlowsThisIndex_[j][cashFlowIndex][0];
562 if (doDeflation_)
563 deflatedCashFlow *= deflatorAndDerivatives_[0];
564 //cashFlowsGenerated_[j][cashFlowIndex].amount[0]*deflatorAndDerivatives_[0];
565 numerairesHeld_[j] += deflatedCashFlow;
566
567 for (Size i=1; i <= numberRates_; ++i)
568 {
569 Real thisDerivative = totalCashFlowsThisIndex_[j][cashFlowIndex][i];
570 if (doDeflation_)
571 {
572 thisDerivative *= deflatorAndDerivatives_[0];
573 thisDerivative += totalCashFlowsThisIndex_[j][cashFlowIndex][0]*deflatorAndDerivatives_[i];
574 fullDerivatives_[i-1] = thisDerivative;
575 }
576 else
577 fullDerivatives_[i-1] = thisDerivative;
578
579 V_[j][stepToUse][i-1] += thisDerivative; // zeroth row of V is t =0 not t_0
580 }
581
582 // ok we've got the derivatives and stored them, now add them to vegas
583 // this corresponds to the \frac{\partial F_n}[\partial theta} term
584 // we add the indirect terms later
585
586 for (Size k=0; k < numberBumps_; ++k)
587 for (Size i=0; i < numberRates_; ++i)
588 {
589 vegasThisPath_[j][k] += fullDerivatives_[i]*jacobiansThisPaths_[stepToUse-1][k][i];
590 }
591
592
593 } // end of (numberCashFlowsThisIndex_[j][cashFlowIndex] > 0)
594 } // end of (Size j=0; j < numberProducts_; ++j)
595 } // end of if (!noFlows)
596 }
597
598 // need to do backwards updating
599 if (flowsFound)
600 {
601 Integer nextStepToUse = std::min<Integer>(currentStep-1, finalStepDone);
602 Integer nextStepIndex = nextStepToUse+1;
603 if (nextStepIndex != stepToUse) // then we need to update V
604 {
605
606 const Matrix& thisPseudoRoot_= pseudoRootStructure_->pseudoRoot(currentStep);
607
608 for (Size i=0; i < numberProducts_; ++i)
609 {
610 // compute partials
611 for (Size f=0; f < factors; ++f)
612 {
613 Real libor = LIBORRates_[stepToUse][numberRates_-1];
614 Real V = V_[i][stepToUse][numberRates_-1];
615 Real pseudo = thisPseudoRoot_[numberRates_-1][f];
616 Real thisPartialTerm = libor*V*pseudo;
617 partials_[f][numberRates_-1] = thisPartialTerm;
618
619 for (Integer r = numberRates_-2; r >=0 ; --r)
620 {
621 Real thisPartialTermr = LIBORRates_[stepToUse][r]*V_[i][stepToUse][r]*thisPseudoRoot_[r][f];
622
623 partials_[f][r] = partials_[f][r+1] + thisPartialTermr;
624
625 }
626 } // end of (Size f=0; f < factors; ++f)
627
628 for (Size j=0; j < numberRates_; ++j)
629 {
630 Real nextV = V_[i][stepToUse][j] * LIBORRatios_[stepToUse][j];
631 V_[i][nextStepIndex][j] = nextV;
632
633 Real summandTerm = 0.0;
634 for (Size f=0; f < factors; ++f)
635 summandTerm += thisPseudoRoot_[j][f]*partials_[f][j];
636
637 summandTerm *= taus[j]*StepsDiscountsSquared_[stepToUse][j];
638
639 V_[i][nextStepIndex][j] += summandTerm;
640
641 } //end of for (Size j=0; j < numberRates_; ++j)
642
643 // we've done the Vs now the vegas
644
645 if (nextStepIndex >0)
646
647 for (Size l=0; l < numberBumps_; ++l)
648 for (Size j=0; j < numberRates_; ++j)
649 vegasThisPath_[i][l] += V_[i][nextStepIndex][j] * jacobiansThisPaths_[nextStepIndex-1][l][j];
650
651
652 } // end of (Size i=0; i < numberProducts_; ++i)
653
654
655
656 } // end of if (nextStepIndex != stepToUse)
657 } // end of if (flowsFound)
658
659 } // end of for (Integer currentStep = numberSteps_-1; currentStep >=0 ; --currentStep)
660
661 // write answer into values
662
663 Size entriesPerProduct = 1+numberRates_+numberBumps_;
664
665 for (Size i=0; i < numberProducts_; ++i)
666 {
667 values[i*entriesPerProduct] = numerairesHeld_[i]*initialNumeraireValue_;
668 for (Size j=0; j < numberRates_; ++j)
669 values[i*entriesPerProduct+1+j] = V_[i][0][j]*initialNumeraireValue_;
670 for (Size k=0; k < numberBumps_; ++k)
671 values[i*entriesPerProduct + numberRates_ +k +1 ] = vegasThisPath_[i][k]*initialNumeraireValue_;
672 }
673
674 return 1.0; // we have put the weight in already, this results in lower variance since weight changes along the path
675 }
676
677 void PathwiseVegasAccountingEngine::multiplePathValues(std::vector<Real>& means, std::vector<Real>& errors,
678 Size numberOfPaths)
679 {
680 std::vector<Real> values(product_->numberOfProducts()*(1+numberRates_+numberBumps_));
681 means.resize(values.size());
682 errors.resize(values.size());
683 std::vector<Real> sums(values.size(),0.0);
684 std::vector<Real> sumsqs(values.size(),0.0);
685
686
687
688 for (Size i=0; i<numberOfPaths; ++i)
689 {
690 /* Real weight = */ singlePathValues(values);
691 // stats.add(values,weight);
692 for (Size j=0; j < values.size(); ++j)
693 {
694 sums[j] += values[j];
695 sumsqs[j] += values[j]*values[j];
696
697 }
698 }
699
700 for (Size j=0; j < values.size(); ++j)
701 {
702 means[j] = sums[j]/numberOfPaths;
703 Real meanSq = sumsqs[j]/numberOfPaths;
704 Real variance = meanSq - means[j]*means[j];
705 errors[j] = std::sqrt(variance/numberOfPaths);
706
707 }
708 }
709
710////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
711////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
712
714 ext::shared_ptr<LogNormalFwdRateEuler> evolver, // method relies heavily on LMM Euler
716 ext::shared_ptr<MarketModel> pseudoRootStructure, // we need pseudo-roots and displacements
717 const std::vector<std::vector<Matrix> >& vegaBumps,
718 Real initialNumeraireValue)
719 : evolver_(std::move(evolver)), product_(product),
720 pseudoRootStructure_(std::move(pseudoRootStructure)), vegaBumps_(vegaBumps),
721 initialNumeraireValue_(initialNumeraireValue), numberProducts_(product->numberOfProducts()),
722 doDeflation_(!product->alreadyDeflated()), numerairesHeld_(product->numberOfProducts()),
723 numberCashFlowsThisStep_(product->numberOfProducts()),
724 cashFlowsGenerated_(product->numberOfProducts()),
725 stepsDiscounts_(pseudoRootStructure_->numberOfRates() + 1),
726 elementary_vegas_ThisPath_(product->numberOfProducts()),
727 deflatorAndDerivatives_(pseudoRootStructure_->numberOfRates() + 1) {
728
729 stepsDiscounts_[0]=1.0;
730
731 numberRates_ = pseudoRootStructure_->numberOfRates();
732 numberSteps_ = pseudoRootStructure_->numberOfSteps();
733 factors_ = pseudoRootStructure_->numberOfFactors();
735
736
737 const EvolutionDescription& evolution = pseudoRootStructure_->evolution();
738 numeraires_ = moneyMarketMeasure(evolution);
739
740
741 QL_REQUIRE(vegaBumps.size() == numberSteps_, "we need precisely one vector of vega bumps for each step.");
742
743 numberBumps_ = vegaBumps[0].size();
744
745 std::vector<Matrix> jacobiansThisPathsModel;
746 for (Size i =0; i < numberRates_; ++i)
747 jacobiansThisPathsModel.emplace_back(numberRates_, factors_);
748
749
750 for (Size i =0; i < numberSteps_; ++i)
751 {
752 jacobianComputers_.emplace_back(
753 pseudoRootStructure_->pseudoRoot(i), evolution.firstAliveRate()[i], numeraires_[i],
754 evolution.rateTaus(), pseudoRootStructure_->displacements());
755
756 // vector of vector of matrices to store jacobians of rates with respect to pseudo-root
757 // elements
758 jacobiansThisPaths_.push_back(jacobiansThisPathsModel);
759 }
760
761
762
764
765
766
768
769 for (Size i=0; i <= numberSteps_; ++i)
770 Discounts_[i][0] = 1.0;
771
772
773 V_.reserve(numberProducts_);
774
775 Matrix modelCashFlowIndex(product_->possibleCashFlowTimes().size(), numberRates_+1);
776
777
779
780 for (Size i=0; i<numberProducts_; ++i)
781 {
782 cashFlowsGenerated_[i].resize(
783 product_->maxNumberOfCashFlowsPerProductPerStep());
784
785 for (auto& j : cashFlowsGenerated_[i])
786 j.amount.resize(numberRates_ + 1);
787
788 numberCashFlowsThisIndex_[i].resize(product_->possibleCashFlowTimes().size());
789
790 V_.push_back(VModel);
791
792
793 totalCashFlowsThisIndex_.push_back(modelCashFlowIndex);
794 }
795
796 LIBORRatios_ = VModel;
797 StepsDiscountsSquared_ = VModel;
798 LIBORRates_ =VModel;
799
800
801
802
803 const std::vector<Time>& cashFlowTimes =
804 product_->possibleCashFlowTimes();
805 numberCashFlowTimes_ = cashFlowTimes.size();
806
807 const std::vector<Time>& rateTimes = product_->evolution().rateTimes();
808 const std::vector<Time>& evolutionTimes = product_->evolution().evolutionTimes();
809 discounters_.reserve(cashFlowTimes.size());
810
811 for (Real cashFlowTime : cashFlowTimes)
812 discounters_.emplace_back(cashFlowTime, rateTimes);
813
814
815 // need to check that we are in money market measure
816
817
818 // we need to allocate cash-flow times to steps, i.e. what is the last step completed before a flow occurs
819 // what we really need is for each step, what cash flow time indices to look at
820
822
823 for (Size i=0; i < numberCashFlowTimes_; ++i)
824 {
825 auto it =
826 std::upper_bound(evolutionTimes.begin(), evolutionTimes.end(), cashFlowTimes[i]);
827 if (it != evolutionTimes.begin())
828 --it;
829 Size index = it - evolutionTimes.begin();
830 cashFlowIndicesThisStep_[index].push_back(i);
831 }
832
834
835// set up this container object
836// std::vector<std::vector<std::vector<Matrix> > > elementary_vegas_ThisPath_; // dimensions are product, step, rate, rate and factor
837
838 { // force destruction of modelVegaMatrix as soon as no longer needed
839 Matrix modelVegaMatrix(numberRates_, factors_,0.0);
840
841 for (Size i=0; i < numberProducts_; ++i)
842 {
844 for (Size j=0; j < numberSteps_; ++j)
845 {
846
847 elementary_vegas_ThisPath_[i][j]= modelVegaMatrix;
848 }
849 }
850 } // modelVegaMatrix destroyed here
851
853/*
854 gaussians_.resize(numberSteps_);
855 distinguishedFactor_=0;
856 distinguishedRate_=0;
857 distinguishedStep_=0;
858*/
859 }
860
862 {
863
864 const std::vector<Real>& initialForwards_(pseudoRootStructure_->initialRates());
865 currentForwards_ = initialForwards_;
866 // clear accumulation variables
867 for (Size i=0; i < numberProducts_; ++i)
868 {
869 numerairesHeld_[i]=0.0;
870
871 for (Size j=0; j < numberCashFlowTimes_; ++j)
872 {
874
875 for (Size k=0; k <= numberRates_; ++k)
876 totalCashFlowsThisIndex_[i][j][k] =0.0;
877 }
878
879 for (Size l=0; l< numberRates_; ++l)
880 for (Size m=0; m <= numberSteps_; ++m)
881 V_[i][m][l] =0.0;
882
883 }
884
885
886
887 Real weight = evolver_->startNewPath();
888 product_->reset();
889
890 Size thisStep;
891
892 bool done = false;
893 do
894 {
895 thisStep = evolver_->currentStep();
896 Size storeStep = thisStep+1;
897 weight *= evolver_->advanceStep();
898
899 done = product_->nextTimeStep(evolver_->currentState(),
902
904 currentForwards_ = evolver_->currentState().forwardRates();
905
906 for (unsigned long i=0; i < numberRates_; ++i)
907 {
908 Real x= evolver_->currentState().discountRatio(i+1,i);
909 stepsDiscounts_[i+1] = x;
910 StepsDiscountsSquared_[storeStep][i] = x*x;
911
912 LIBORRatios_[storeStep][i] = currentForwards_[i]/lastForwards_[i];
913 LIBORRates_[storeStep][i] = currentForwards_[i];
914 Discounts_[storeStep][i+1] = evolver_->currentState().discountRatio(i+1,0);
915 }
916
917 jacobianComputers_[thisStep].getBumps(lastForwards_,
920 evolver_->browniansThisStep(),
921 jacobiansThisPaths_[thisStep]);
922
923// gaussians_[thisStep] = evolver_->browniansThisStep();
924
925
926
927 // for each product...
928 for (Size i=0; i<numberProducts_; ++i)
929 {
930 // ...and each cash flow...
931 for (Size j=0; j<numberCashFlowsThisStep_[i]; ++j)
932 {
933 Size k = cashFlowsGenerated_[i][j].timeIndex;
935
936 for (Size l=0; l <= numberRates_; ++l)
937 totalCashFlowsThisIndex_[i][k][l] += cashFlowsGenerated_[i][j].amount[l]*weight;
938
939 }
940 }
941
942
943 } while (!done);
944
945 // ok we've gathered cash-flows, still have to backwards computation
946
947 Size factors = pseudoRootStructure_->numberOfFactors();
948 const std::vector<Time>& taus= pseudoRootStructure_->evolution(). rateTaus();
949
950 bool flowsFound = false;
951
952 Integer finalStepDone = thisStep;
953
954 for (Integer currentStep = numberSteps_-1; currentStep >=0 ; --currentStep) // must be a signed type as we go negative
955 {
956 Integer stepToUse = std::min<Integer>(currentStep, finalStepDone)+1;
957
958 for (Size k=0; k < cashFlowIndicesThisStep_[currentStep].size(); ++k)
959 {
960 Size cashFlowIndex =cashFlowIndicesThisStep_[currentStep][k];
961
962 // first check to see if anything actually happened before spending time on computing stuff
963 bool noFlows = true;
964 for (Size l=0; l < numberProducts_ && noFlows; ++l)
965 noFlows = noFlows && (numberCashFlowsThisIndex_[l][cashFlowIndex] ==0);
966
967 flowsFound = flowsFound || !noFlows;
968
969 if (!noFlows)
970 {
971 if (doDeflation_)
972 discounters_[cashFlowIndex].getFactors(LIBORRates_, Discounts_,stepToUse, deflatorAndDerivatives_); // get amount to discount cash flow by and amount to multiply its derivatives by
973
974 for (Size j=0; j < numberProducts_; ++j)
975 {
976 if (numberCashFlowsThisIndex_[j][cashFlowIndex] > 0)
977 {
978 Real deflatedCashFlow = totalCashFlowsThisIndex_[j][cashFlowIndex][0];
979 if (doDeflation_)
980 deflatedCashFlow *= deflatorAndDerivatives_[0];
981 //cashFlowsGenerated_[j][cashFlowIndex].amount[0]*deflatorAndDerivatives_[0];
982 numerairesHeld_[j] += deflatedCashFlow;
983
984 for (Size i=1; i <= numberRates_; ++i)
985 {
986 Real thisDerivative = totalCashFlowsThisIndex_[j][cashFlowIndex][i];
987 if (doDeflation_)
988 {
989 thisDerivative *= deflatorAndDerivatives_[0];
990 thisDerivative += totalCashFlowsThisIndex_[j][cashFlowIndex][0]*deflatorAndDerivatives_[i];
991 fullDerivatives_[i-1] = thisDerivative;
992 }
993 else
994 fullDerivatives_[i-1] = thisDerivative;
995
996 V_[j][stepToUse][i-1] += thisDerivative; // zeroth row of V is t =0 not t_0
997 } // end of for (Size i=1; i <= numberRates_; ++i)
998 } // end of (numberCashFlowsThisIndex_[j][cashFlowIndex] > 0)
999 } // end of (Size j=0; j < numberProducts_; ++j)
1000 } // end of if (!noFlows)
1001 }
1002
1003 // need to do backwards updating
1004 if (flowsFound)
1005 {
1006 Integer nextStepToUse = std::min<Integer>(currentStep-1, finalStepDone);
1007 Integer nextStepIndex = nextStepToUse+1;
1008 if (nextStepIndex != stepToUse) // then we need to update V
1009 {
1010
1011 const Matrix& thisPseudoRoot_= pseudoRootStructure_->pseudoRoot(currentStep);
1012
1013 for (Size i=0; i < numberProducts_; ++i)
1014 {
1015 // compute partials
1016 for (Size f=0; f < factors; ++f)
1017 {
1018 Real libor = LIBORRates_[stepToUse][numberRates_-1];
1019 Real V = V_[i][stepToUse][numberRates_-1];
1020 Real pseudo = thisPseudoRoot_[numberRates_-1][f];
1021 Real thisPartialTerm = libor*V*pseudo;
1022 partials_[f][numberRates_-1] = thisPartialTerm;
1023
1024 for (Integer r = numberRates_-2; r >=0 ; --r)
1025 {
1026 Real thisPartialTermr = LIBORRates_[stepToUse][r]*V_[i][stepToUse][r]*thisPseudoRoot_[r][f];
1027
1028 partials_[f][r] = partials_[f][r+1] + thisPartialTermr;
1029
1030 }
1031 } // end of (Size f=0; f < factors; ++f)
1032
1033 for (Size j=0; j < numberRates_; ++j)
1034 {
1035 Real nextV = V_[i][stepToUse][j] * LIBORRatios_[stepToUse][j];
1036 V_[i][nextStepIndex][j] = nextV;
1037
1038 Real summandTerm = 0.0;
1039 for (Size f=0; f < factors; ++f)
1040 summandTerm += thisPseudoRoot_[j][f]*partials_[f][j];
1041
1042 summandTerm *= taus[j]*StepsDiscountsSquared_[stepToUse][j];
1043
1044 V_[i][nextStepIndex][j] += summandTerm;
1045
1046 } //end of for (Size j=0; j < numberRates_; ++j)
1047
1048 } // end of (Size i=0; i < numberProducts_; ++i)
1049
1050 } // end of if (nextStepIndex != stepToUse)
1051
1052 } // end of if (flowsFound)
1053
1054 } // end of for (Integer currentStep = numberSteps_-1; currentStep >=0 ; --currentStep)
1055
1056
1057 // all V matrices computed we now compute the elementary vegas for this path
1058
1059 for (Size i=0; i < numberProducts_; ++i)
1060 {
1061 for (Size j=0; j < numberSteps_; ++j)
1062 {
1063 Size nextIndex = j+1;
1064
1065 // we know V, we need to pair against the senstivity of the rate to the elementary vega
1066 // note the simplification here arising from the fact that the elementary vega affects the evolution on precisely one step
1067
1068 for (Size k=0; k < numberRates_; ++k)
1069 for (Size f=0; f < factors_; ++f)
1070 {
1071 Real sensitivity =0.0;
1072
1073 for (Size r=0; r < numberRates_; ++r)
1074 {
1075 sensitivity += V_[i][nextIndex][r]*jacobiansThisPaths_[j][r][k][f];
1076
1077 }
1078/*
1079 if (j ==distinguishedStep_ && k ==distinguishedRate_ &&f== distinguishedFactor_)
1080 std::cout << sensitivity << "," << jacobiansThisPaths_[j][j][k][f] << "," << gaussians_[j][f] << "," << V_[i][nextIndex][j] << "," << LIBORRates_[nextIndex][j] << "\n";
1081 */
1082
1083
1084 elementary_vegas_ThisPath_[i][j][k][f] = sensitivity;
1085 }
1086 }
1087 }
1088
1089
1090 // write answer into values
1091
1092 Size entriesPerProduct = 1+numberRates_+numberElementaryVegas_;
1093
1094 for (Size i=0; i < numberProducts_; ++i)
1095 {
1096 values[i*entriesPerProduct] = numerairesHeld_[i]*initialNumeraireValue_;
1097 for (Size j=0; j < numberRates_; ++j)
1098 values[i*entriesPerProduct+1+j] = V_[i][0][j]*initialNumeraireValue_;
1099
1100 for (Size k=0; k < numberSteps_; ++k)
1101 for (Size l=0; l < numberRates_; ++l)
1102 for (Size m=0; m < factors_; ++m)
1103 values[i*entriesPerProduct + numberRates_ +1 + m+ l*factors_ + k*numberRates_*factors_] = elementary_vegas_ThisPath_[i][k][l][m]*initialNumeraireValue_;
1104
1105 }
1106
1107 return 1.0; // we have put the weight in already, this results in lower variance since weight changes along the path
1108
1109}
1110
1111 void PathwiseVegasOuterAccountingEngine::multiplePathValuesElementary(std::vector<Real>& means, std::vector<Real>& errors,
1112 Size numberOfPaths)
1113 {
1114 Size numberOfElementaryVegas = numberRates_*numberSteps_*factors_;
1115
1116 std::vector<Real> values(product_->numberOfProducts()*(1+numberRates_+numberOfElementaryVegas));
1117 means.resize(values.size());
1118 errors.resize(values.size());
1119 std::vector<Real> sums(values.size(),0.0);
1120 std::vector<Real> sumsqs(values.size(),0.0);
1121
1122
1123
1124 for (Size i=0; i<numberOfPaths; ++i)
1125 {
1126 singlePathValues(values);
1127
1128 for (Size j=0; j < values.size(); ++j)
1129 {
1130 sums[j] += values[j];
1131 sumsqs[j] += values[j]*values[j];
1132
1133 }
1134 }
1135
1136 for (Size j=0; j < values.size(); ++j)
1137 {
1138 means[j] = sums[j]/numberOfPaths;
1139 Real meanSq = sumsqs[j]/numberOfPaths;
1140 Real variance = meanSq - means[j]*means[j];
1141 errors[j] = std::sqrt(variance/numberOfPaths);
1142
1143 }
1144 }
1145
1146 void PathwiseVegasOuterAccountingEngine::multiplePathValues(std::vector<Real>& means, std::vector<Real>& errors,Size numberOfPaths)
1147 {
1148 std::vector<Real> allMeans;
1149 std::vector<Real> allErrors;
1150
1151 multiplePathValuesElementary(allMeans,allErrors,numberOfPaths);
1152
1153 Size outDataPerProduct = 1+numberRates_+numberBumps_;
1154 Size inDataPerProduct = 1+numberRates_+numberElementaryVegas_;
1155
1156 means.resize((1+numberRates_+numberBumps_)*numberProducts_);
1157 errors.resize((1+numberRates_+numberBumps_)*numberProducts_); // post linear combinations, errors are not meaningful so don't attempt to compute s.e.s for vegas
1158
1159 for (Size p=0; p < numberProducts_; ++p)
1160 {
1161 for (Size i=0; i < 1 + numberRates_; ++i)
1162 {
1163 means[i+p*outDataPerProduct] = allMeans[i+p*inDataPerProduct];
1164 errors[i+p*outDataPerProduct] = allErrors[i+p*inDataPerProduct];
1165 }
1166
1167 for (Size bump=0; bump<numberBumps_; ++bump)
1168 {
1169 Real thisVega=0.0;
1170
1171
1172 for (Size t=0; t < numberSteps_; ++t)
1173 for (Size r=0; r < numberRates_; ++r)
1174 for (Size f=0; f < factors_; ++f)
1175 thisVega+= vegaBumps_[t][bump][r][f]*allMeans[p*inDataPerProduct+1+numberRates_+t*numberRates_*factors_+r*factors_+f];
1176
1177
1178 means[p*outDataPerProduct+1+numberRates_+bump] = thisVega;
1179 }
1180
1181 }
1182
1183 } // end of method
1184
1185} // end of namespace
1186
1187
1188
1189
1190
1191
1192
1193
1194
cloning proxy to an underlying object
Definition: clone.hpp:40
Market-model evolution description.
const std::vector< Time > & rateTaus() const
const std::vector< Size > & firstAliveRate() const
Statistics analysis of N-dimensional (sequence) data.
void add(const Sequence &sample, Real weight=1.0)
Matrix used in linear algebra.
Definition: matrix.hpp:41
ext::shared_ptr< LogNormalFwdRateEuler > evolver_
std::vector< std::vector< MarketModelPathwiseMultiProduct::CashFlow > > cashFlowsGenerated_
PathwiseAccountingEngine(ext::shared_ptr< LogNormalFwdRateEuler > evolver, const Clone< MarketModelPathwiseMultiProduct > &product, ext::shared_ptr< MarketModel > pseudoRootStructure, Real initialNumeraireValue)
std::vector< std::vector< Size > > numberCashFlowsThisIndex_
void multiplePathValues(SequenceStatisticsInc &stats, Size numberOfPaths)
ext::shared_ptr< MarketModel > pseudoRootStructure_
std::vector< std::vector< Size > > cashFlowIndicesThisStep_
Real singlePathValues(std::vector< Real > &values)
Clone< MarketModelPathwiseMultiProduct > product_
std::vector< MarketModelPathwiseDiscounter > discounters_
void multiplePathValues(std::vector< Real > &means, std::vector< Real > &errors, Size numberOfPaths)
PathwiseVegasAccountingEngine(ext::shared_ptr< LogNormalFwdRateEuler > evolver, const Clone< MarketModelPathwiseMultiProduct > &product, ext::shared_ptr< MarketModel > pseudoRootStructure, const std::vector< std::vector< Matrix > > &VegaBumps, Real initialNumeraireValue)
ext::shared_ptr< LogNormalFwdRateEuler > evolver_
std::vector< std::vector< MarketModelPathwiseMultiProduct::CashFlow > > cashFlowsGenerated_
std::vector< std::vector< Size > > numberCashFlowsThisIndex_
std::vector< RatePseudoRootJacobian > jacobianComputers_
ext::shared_ptr< MarketModel > pseudoRootStructure_
std::vector< std::vector< Size > > cashFlowIndicesThisStep_
Real singlePathValues(std::vector< Real > &values)
Clone< MarketModelPathwiseMultiProduct > product_
std::vector< MarketModelPathwiseDiscounter > discounters_
void multiplePathValues(std::vector< Real > &means, std::vector< Real > &errors, Size numberOfPaths)
Use to get vegas with respect to VegaBumps.
ext::shared_ptr< LogNormalFwdRateEuler > evolver_
void multiplePathValuesElementary(std::vector< Real > &means, std::vector< Real > &errors, Size numberOfPaths)
Use to get vegas with respect to pseudo-root-elements.
std::vector< std::vector< MarketModelPathwiseMultiProduct::CashFlow > > cashFlowsGenerated_
std::vector< std::vector< Size > > numberCashFlowsThisIndex_
std::vector< std::vector< Matrix > > vegaBumps_
std::vector< std::vector< Matrix > > jacobiansThisPaths_
std::vector< std::vector< Matrix > > elementary_vegas_ThisPath_
PathwiseVegasOuterAccountingEngine(ext::shared_ptr< LogNormalFwdRateEuler > evolver, const Clone< MarketModelPathwiseMultiProduct > &product, ext::shared_ptr< MarketModel > pseudoRootStructure, const std::vector< std::vector< Matrix > > &VegaBumps, Real initialNumeraireValue)
std::vector< RatePseudoRootJacobianAllElements > jacobianComputers_
std::vector< std::vector< Size > > cashFlowIndicesThisStep_
Clone< MarketModelPathwiseMultiProduct > product_
std::vector< MarketModelPathwiseDiscounter > discounters_
const DefaultType & t
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
Definition: errors.hpp:117
LinearInterpolation variance
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
std::vector< Size > moneyMarketMeasure(const EvolutionDescription &evol)
STL namespace.
ext::shared_ptr< YieldTermStructure > r
std::vector< Size > numberCashFlowsThisStep_
std::vector< std::vector< CashFlow > > cashFlowsGenerated_