1 /* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ 2 3 /* 4 Copyright (C) 2006 Ferdinando Ametrano 5 Copyright (C) 2006 Mark Joshi 6 7 This file is part of QuantLib, a free-software/open-source library 8 for financial quantitative analysts and developers - http://quantlib.org/ 9 10 QuantLib is free software: you can redistribute it and/or modify it 11 under the terms of the QuantLib license. You should have received a 12 copy 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 16 This program is distributed in the hope that it will be useful, but WITHOUT 17 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 18 FOR A PARTICULAR PURPOSE. See the license for more details. 19 */ 20 21 #include <ql/models/marketmodels/evolvers/lognormalfwdrateipc.hpp> 22 #include <ql/models/marketmodels/marketmodel.hpp> 23 #include <ql/models/marketmodels/evolutiondescription.hpp> 24 #include <ql/models/marketmodels/browniangenerator.hpp> 25 #include <ql/models/marketmodels/driftcomputation/lmmdriftcalculator.hpp> 26 27 namespace QuantLib { 28 LogNormalFwdRateIpc(const ext::shared_ptr<MarketModel> & marketModel,const BrownianGeneratorFactory & factory,const std::vector<Size> & numeraires,Size initialStep)29 LogNormalFwdRateIpc::LogNormalFwdRateIpc( 30 const ext::shared_ptr<MarketModel>& marketModel, 31 const BrownianGeneratorFactory& factory, 32 const std::vector<Size>& numeraires, 33 Size initialStep) 34 : marketModel_(marketModel), 35 numeraires_(numeraires), 36 initialStep_(initialStep), 37 numberOfRates_(marketModel->numberOfRates()), 38 numberOfFactors_(marketModel->numberOfFactors()), 39 curveState_(marketModel->evolution().rateTimes()), 40 forwards_(marketModel->initialRates()), 41 displacements_(marketModel->displacements()), 42 logForwards_(numberOfRates_), initialLogForwards_(numberOfRates_), 43 drifts1_(numberOfRates_), initialDrifts_(numberOfRates_), 44 g_(numberOfRates_), brownians_(numberOfFactors_), 45 correlatedBrownians_(numberOfRates_), 46 rateTaus_(marketModel->evolution().rateTaus()), 47 alive_(marketModel->evolution().firstAliveRate()) 48 { 49 checkCompatibility(marketModel->evolution(), numeraires); 50 QL_REQUIRE(isInTerminalMeasure(marketModel->evolution(), numeraires), 51 "terminal measure required for ipc "); 52 53 Size steps = marketModel->evolution().numberOfSteps(); 54 55 generator_ = factory.create(numberOfFactors_, steps-initialStep_); 56 57 currentStep_ = initialStep_; 58 59 calculators_.reserve(steps); 60 fixedDrifts_.reserve(steps); 61 for (Size j=0; j<steps; ++j) { 62 const Matrix& A = marketModel->pseudoRoot(j); 63 calculators_.push_back(LMMDriftCalculator(A, 64 displacements_, 65 marketModel->evolution().rateTaus(), 66 numeraires[j], 67 alive_[j])); 68 const Matrix& C = marketModel->covariance(j); 69 std::vector<Real> fixed(numberOfRates_); 70 for (Size k=0; k<numberOfRates_; ++k) { 71 Real variance = C[k][k]; 72 fixed[k] = -0.5*variance; 73 } 74 fixedDrifts_.push_back(fixed); 75 } 76 77 setForwards(marketModel_->initialRates()); 78 } 79 numeraires() const80 const std::vector<Size>& LogNormalFwdRateIpc::numeraires() const { 81 return numeraires_; 82 } 83 setForwards(const std::vector<Real> & forwards)84 void LogNormalFwdRateIpc::setForwards(const std::vector<Real>& forwards) 85 { 86 QL_REQUIRE(forwards.size()==numberOfRates_, 87 "mismatch between forwards and rateTimes"); 88 for (Size i=0; i<numberOfRates_; ++i) 89 initialLogForwards_[i] = std::log(forwards[i] + 90 displacements_[i]); 91 calculators_[initialStep_].compute(forwards, initialDrifts_); 92 } 93 setInitialState(const CurveState & cs)94 void LogNormalFwdRateIpc::setInitialState(const CurveState& cs) { 95 setForwards(cs.forwardRates()); 96 } 97 startNewPath()98 Real LogNormalFwdRateIpc::startNewPath() { 99 currentStep_ = initialStep_; 100 std::copy(initialLogForwards_.begin(), initialLogForwards_.end(), 101 logForwards_.begin()); 102 return generator_->nextPath(); 103 } 104 advanceStep()105 Real LogNormalFwdRateIpc::advanceStep() 106 { 107 // we're going from T1 to T2: 108 109 // a) compute drifts D1 at T1; 110 if (currentStep_ > initialStep_) { 111 calculators_[currentStep_].computePlain(forwards_, drifts1_); 112 } else { 113 std::copy(initialDrifts_.begin(), initialDrifts_.end(), 114 drifts1_.begin()); 115 } 116 117 Real weight = generator_->nextStep(brownians_); 118 const Matrix& A = marketModel_->pseudoRoot(currentStep_); 119 const Matrix& C = marketModel_->covariance(currentStep_); 120 const std::vector<Real>& fixedDrift = fixedDrifts_[currentStep_]; 121 122 Integer alive = alive_[currentStep_]; 123 Real drifts2; 124 for (Integer i=numberOfRates_-1; i>=alive; --i) { 125 drifts2 = 0.0; 126 for (Size j=i+1; j<numberOfRates_; ++j) { 127 drifts2 -= g_[j]*C[i][j]; 128 } 129 logForwards_[i] += 0.5*(drifts1_[i]+drifts2) + fixedDrift[i]; 130 logForwards_[i] += 131 std::inner_product(A.row_begin(i), A.row_end(i), 132 brownians_.begin(), 0.0); 133 forwards_[i] = std::exp(logForwards_[i]) - displacements_[i]; 134 g_[i] = rateTaus_[i]*(forwards_[i]+displacements_[i])/ 135 (1.0+rateTaus_[i]*forwards_[i]); 136 } 137 138 // update curve state 139 curveState_.setOnForwardRates(forwards_); 140 141 ++currentStep_; 142 143 return weight; 144 } 145 currentStep() const146 Size LogNormalFwdRateIpc::currentStep() const { 147 return currentStep_; 148 } 149 currentState() const150 const CurveState& LogNormalFwdRateIpc::currentState() const { 151 return curveState_; 152 } 153 154 } 155