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