1 /* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ 2 3 /* 4 Copyright (C) 2020 Lew Wei Hao 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/methods/finitedifferences/meshers/fdmmesher.hpp> 21 #include <ql/methods/finitedifferences/operators/fdmcirop.hpp> 22 #include <ql/methods/finitedifferences/operators/fdmlinearoplayout.hpp> 23 #include <ql/methods/finitedifferences/operators/secondderivativeop.hpp> 24 #include <ql/methods/finitedifferences/operators/secondordermixedderivativeop.hpp> 25 #include <ql/processes/blackscholesprocess.hpp> 26 27 namespace QuantLib { 28 FdmCIREquityPart(const ext::shared_ptr<FdmMesher> & mesher,const ext::shared_ptr<GeneralizedBlackScholesProcess> & bsProcess,Real strike)29 FdmCIREquityPart::FdmCIREquityPart( 30 const ext::shared_ptr<FdmMesher>& mesher, 31 const ext::shared_ptr<GeneralizedBlackScholesProcess>& bsProcess, 32 Real strike) 33 : dxMap_ (FirstDerivativeOp(0, mesher)), 34 dxxMap_(SecondDerivativeOp(0, mesher)), 35 mapT_ (0, mesher), 36 mesher_(mesher), 37 qTS_(bsProcess->dividendYield().currentLink()), 38 strike_(strike), 39 sigma1_(bsProcess->blackVolatility().currentLink()){ 40 } 41 setTime(Time t1,Time t2)42 void FdmCIREquityPart::setTime(Time t1, Time t2) { 43 const Rate q = qTS_->forwardRate(t1, t2, Continuous).rate(); 44 45 const Real v = sigma1_->blackForwardVariance(t1, t2, strike_)/(t2-t1); 46 47 mapT_.axpyb(mesher_->locations(1) - q - 0.5*v, dxMap_, 48 dxxMap_.mult(Array(mesher_->layout()->size(), v/2)), -0.5*mesher_->locations(1)); 49 } 50 getMap() const51 const TripleBandLinearOp& FdmCIREquityPart::getMap() const { 52 return mapT_; 53 } 54 FdmCIRRatesPart(const ext::shared_ptr<FdmMesher> & mesher,Real sigma,Real kappa,Real theta)55 FdmCIRRatesPart::FdmCIRRatesPart( 56 const ext::shared_ptr<FdmMesher>& mesher, 57 Real sigma, Real kappa, Real theta) 58 : dyMap_(SecondDerivativeOp(1, mesher) 59 .mult(sigma*sigma*mesher->locations(1)) 60 .add(FirstDerivativeOp(1, mesher) 61 .mult(kappa*(theta - mesher->locations(1))))), 62 mapT_(1, mesher), 63 mesher_(mesher){ 64 } 65 setTime(Time t1,Time t2)66 void FdmCIRRatesPart::setTime(Time t1, Time t2) { 67 mapT_.axpyb(Array(), dyMap_, dyMap_, -0.5*mesher_->locations(1)); 68 } 69 getMap() const70 const TripleBandLinearOp& FdmCIRRatesPart::getMap() const { 71 return mapT_; 72 } 73 FdmCIRMixedPart(const ext::shared_ptr<FdmMesher> & mesher,const ext::shared_ptr<CoxIngersollRossProcess> & cirProcess,const ext::shared_ptr<GeneralizedBlackScholesProcess> & bsProcess,const Real rho,const Real strike)74 FdmCIRMixedPart::FdmCIRMixedPart( 75 const ext::shared_ptr<FdmMesher>& mesher, 76 const ext::shared_ptr<CoxIngersollRossProcess> & cirProcess, 77 const ext::shared_ptr<GeneralizedBlackScholesProcess> & bsProcess, 78 const Real rho, 79 const Real strike) 80 : dyMap_(SecondOrderMixedDerivativeOp(0, 1, mesher) 81 .mult(Array(mesher->layout()->size(), 2*rho*cirProcess->volatility()))), 82 mapT_(0, 1, mesher), 83 mesher_(mesher), 84 sigma1_(bsProcess->blackVolatility().currentLink()), 85 strike_(strike){ 86 } 87 setTime(Time t1,Time t2)88 void FdmCIRMixedPart::setTime(Time t1, Time t2) { 89 const Real v = std::sqrt(sigma1_->blackForwardVariance(t1, t2, strike_)/(t2-t1)); 90 NinePointLinearOp op(dyMap_.mult(Array(mesher_->layout()->size(), v))); 91 mapT_.swap(op); 92 } 93 getMap() const94 const NinePointLinearOp& FdmCIRMixedPart::getMap() const { 95 return mapT_; 96 } 97 FdmCIROp(const ext::shared_ptr<FdmMesher> & mesher,const ext::shared_ptr<CoxIngersollRossProcess> & cirProcess,const ext::shared_ptr<GeneralizedBlackScholesProcess> & bsProcess,const Real rho,const Real strike)98 FdmCIROp::FdmCIROp( 99 const ext::shared_ptr<FdmMesher>& mesher, 100 const ext::shared_ptr<CoxIngersollRossProcess> & cirProcess, 101 const ext::shared_ptr<GeneralizedBlackScholesProcess> & bsProcess, 102 const Real rho, 103 const Real strike) 104 : dxMap_(mesher, 105 bsProcess, 106 strike), 107 dyMap_(mesher, 108 cirProcess->volatility(), 109 cirProcess->speed(), 110 cirProcess->level()), 111 dzMap_(mesher, 112 cirProcess, 113 bsProcess, 114 rho, 115 strike){ 116 } 117 118 setTime(Time t1,Time t2)119 void FdmCIROp::setTime(Time t1, Time t2) { 120 dxMap_.setTime(t1, t2); 121 dyMap_.setTime(t1, t2); 122 dzMap_.setTime(t1, t2); 123 } 124 size() const125 Size FdmCIROp::size() const { 126 return 2; 127 } 128 apply(const Array & u) const129 Disposable<Array> FdmCIROp::apply(const Array& u) const { 130 Array dx = dxMap_.getMap().apply(u); 131 Array dy = dyMap_.getMap().apply(u); 132 Array dz = dzMap_.getMap().apply(u); 133 134 return (dy + dx + dz); 135 } 136 apply_direction(Size direction,const Array & r) const137 Disposable<Array> FdmCIROp::apply_direction(Size direction, 138 const Array& r) const { 139 if (direction == 0) 140 return dxMap_.getMap().apply(r); 141 else if (direction == 1) 142 return dyMap_.getMap().apply(r); 143 else 144 QL_FAIL("direction too large"); 145 } 146 apply_mixed(const Array & r) const147 Disposable<Array> FdmCIROp::apply_mixed(const Array& r) const { 148 return dzMap_.getMap().apply(r); 149 } 150 151 Disposable<Array> solve_splitting(Size direction,const Array & r,Real a) const152 FdmCIROp::solve_splitting(Size direction, 153 const Array& r, Real a) const { 154 155 if (direction == 0) { 156 return dxMap_.getMap().solve_splitting(r, a, 1.0); 157 } 158 else if (direction == 1) { 159 return dyMap_.getMap().solve_splitting(r, a, 1.0); 160 } 161 else 162 QL_FAIL("direction too large"); 163 } 164 165 Disposable<Array> preconditioner(const Array & r,Real dt) const166 FdmCIROp::preconditioner(const Array& r, Real dt) const { 167 168 return solve_splitting(1, solve_splitting(0, r, dt), dt) ; 169 } 170 171 #if !defined(QL_NO_UBLAS_SUPPORT) 172 Disposable<std::vector<SparseMatrix> > toMatrixDecomp() const173 FdmCIROp::toMatrixDecomp() const { 174 std::vector<SparseMatrix> retVal(3); 175 176 retVal[0] = dxMap_.getMap().toMatrix(); 177 retVal[1] = dyMap_.getMap().toMatrix(); 178 retVal[2] = dzMap_.getMap().toMatrix(); 179 180 return retVal; 181 } 182 #endif 183 } 184