QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
fdmhestonfwdop.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) 2012, 2013 Klaus Spanderen
5 Copyright (C) 2014 Johannes Göttker-Schnetmann
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/*! \file fdmhestonfwdop.cpp
22*/
23
32#include <cmath>
33#include <utility>
34
35using std::exp;
36
37namespace QuantLib {
38
39 FdmHestonFwdOp::FdmHestonFwdOp(const ext::shared_ptr<FdmMesher>& mesher,
40 const ext::shared_ptr<HestonProcess>& process,
42 ext::shared_ptr<LocalVolTermStructure> leverageFct,
43 const Real mixingFactor)
44 : type_(type), kappa_(process->kappa()), theta_(process->theta()), sigma_(process->sigma()),
45 rho_(process->rho()), v0_(process->v0()), mixedSigma_(mixingFactor * sigma_),
46 rTS_(process->riskFreeRate().currentLink()), qTS_(process->dividendYield().currentLink()),
47 varianceValues_(0.5 * mesher->locations(1)),
48 dxMap_(ext::make_shared<FirstDerivativeOp>(0, mesher)),
49 dxxMap_(ext::make_shared<ModTripleBandLinearOp>(TripleBandLinearOp(
50 type == FdmSquareRootFwdOp::Log ?
51 SecondDerivativeOp(0, mesher).mult(0.5 * Exp(mesher->locations(1))) :
52 SecondDerivativeOp(0, mesher).mult(0.5 * mesher->locations(1))))),
53 boundary_(ext::make_shared<ModTripleBandLinearOp>(TripleBandLinearOp(
54 SecondDerivativeOp(0, mesher).mult(Array(mesher->locations(0).size(), 0.0))))),
55 mapX_(ext::make_shared<TripleBandLinearOp>(0, mesher)),
56 mapY_(ext::make_shared<FdmSquareRootFwdOp>(mesher, kappa_, theta_, mixedSigma_, 1, type)),
57 correlation_(ext::make_shared<NinePointLinearOp>(
58 type == FdmSquareRootFwdOp::Log ?
60 .mult(Array(mesher->layout()->size(), rho_ * mixedSigma_)) :
62 .mult(rho_ * mixedSigma_ * mesher->locations(1)))),
63 leverageFct_(std::move(leverageFct)), mesher_(mesher) {
64 // zero flux boundary condition
65 const Size n = mesher->layout()->dim()[1];
66 const Real lowerBoundaryFactor = mapY_->lowerBoundaryFactor(type);
67 const Real upperBoundaryFactor = mapY_->upperBoundaryFactor(type);
68
69 const Real logFacLow = type == FdmSquareRootFwdOp::Log ? Real(exp(mapY_->v(0))) : 1.0;
70 const Real logFacUpp = type == FdmSquareRootFwdOp::Log ? Real(exp(mapY_->v(n+1))) : 1.0;
71
72 const Real alpha = -2*rho_/mixedSigma_*lowerBoundaryFactor*logFacLow;
73 const Real beta = -2*rho_/mixedSigma_*upperBoundaryFactor*logFacUpp;
74
76
77 for (const auto& iter : *mesher->layout()) {
78 if (iter.coordinates()[1] == 0) {
79 const Size idx = iter.index();
80 if (!leverageFct_) {
81 dxxMap_->upper(idx) += alpha*fDx.upper(idx);
82 dxxMap_->diag(idx) += alpha*fDx.diag(idx);
83 dxxMap_->lower(idx) += alpha*fDx.lower(idx);
84 }
85 boundary_->upper(idx)= alpha*fDx.upper(idx);
86 boundary_->diag(idx) = alpha*fDx.diag(idx);
87 boundary_->lower(idx) = alpha*fDx.lower(idx);
88 }
89 else if (iter.coordinates()[1] == n-1) {
90 const Size idx = iter.index();
91
92 if (!leverageFct_) {
93 dxxMap_->upper(idx)+= beta*fDx.upper(idx);
94 dxxMap_->diag(idx) += beta*fDx.diag(idx);
95 dxxMap_->lower(idx) += beta*fDx.lower(idx);
96 }
97 boundary_->upper(idx)= beta*fDx.upper(idx);
98 boundary_->diag(idx) = beta*fDx.diag(idx);
99 boundary_->lower(idx) = beta*fDx.lower(idx);
100 }
101 }
102 }
103
105 return 2;
106 }
107
109 const Rate r = rTS_->forwardRate(t1, t2, Continuous).rate();
110 const Rate q = qTS_->forwardRate(t1, t2, Continuous).rate();
111 if (leverageFct_ != nullptr) {
112 L_ = getLeverageFctSlice(t1, t2);
113 Array Lsquare = L_*L_;
115 mapX_->axpyb( Array(1, -r + q), *dxMap_,
116 dxxMap_->multR(Lsquare).add(boundary_->multR(L_))
117 .add(dxMap_->multR(rho_*mixedSigma_*L_))
118 .add(dxMap_->mult(varianceValues_).multR(Lsquare)),
119 Array());
120 } else if (type_ == FdmSquareRootFwdOp::Power) {
121 mapX_->axpyb( Array(1, -r + q), *dxMap_,
122 dxxMap_->multR(Lsquare).add(boundary_->multR(L_))
123 .add(dxMap_->multR(rho_*2.0*kappa_*theta_/(mixedSigma_)*L_))
124 .add(dxMap_->mult(varianceValues_).multR(Lsquare)), Array());
125 } else if (type_ == FdmSquareRootFwdOp::Log) {
126 mapX_->axpyb( Array(1, -r + q), *dxMap_,
127 dxxMap_->multR(Lsquare).add(boundary_->multR(L_))
128 .add(dxMap_->mult(0.5*Exp(2.0*varianceValues_)).multR(Lsquare)),
129 Array());
130 }
131 } else {
133 mapX_->axpyb( - r + q + rho_*mixedSigma_ + varianceValues_, *dxMap_,
134 *dxxMap_, Array());
135 } else if (type_ == FdmSquareRootFwdOp::Power) {
136 mapX_->axpyb( - r + q + rho_*2.0*kappa_*theta_/(mixedSigma_) + varianceValues_,
137 *dxMap_, *dxxMap_, Array());
138 } else if (type_ == FdmSquareRootFwdOp::Log) {
139 mapX_->axpyb( - r + q + 0.5*Exp(2.0*varianceValues_), *dxMap_,
140 *dxxMap_, Array());
141 }
142 }
143 }
144
146 if (leverageFct_ != nullptr) {
147 return mapX_->apply(u)
148 + mapY_->apply(u)
149 + correlation_->apply(L_*u);
150 } else {
151 return mapX_->apply(u)
152 + mapY_->apply(u)
153 + correlation_->apply(u);
154 }
155 }
156
158 if (leverageFct_ != nullptr) {
159 return correlation_->apply(L_*u);
160 } else {
161 return correlation_->apply(u);
162 }
163 }
164
166 Size direction, const Array& u) const {
167
168 if (direction == 0)
169 return mapX_->apply(u) ;
170 else if (direction == 1)
171 return mapY_->apply(u) ;
172 else
173 QL_FAIL("direction too large");
174 }
175
177 Size direction, const Array& u, Real s) const{
178 if (direction == 0) {
179 return mapX_->solve_splitting(u, s, 1.0);
180 }
181 else if (direction == 1) {
182 return mapY_->solve_splitting(1, u, s);
183 }
184 else
185 QL_FAIL("direction too large");
186 }
187
189 const Array& u, Real dt) const{
190 return solve_splitting(1, u, dt);
191 }
192
194 Array v(mesher_->layout()->size(), 1.0);
195
196 if (!leverageFct_)
197 return v;
198
199 const Real t = 0.5*(t1+t2);
200 const Time time = std::min(leverageFct_->maxTime(), t);
201 //std::max(leverageFct_->minTime(), t));
202
203 for (const auto& iter : *mesher_->layout()) {
204 const Size nx = iter.coordinates()[0];
205
206 if (iter.coordinates()[1] == 0) {
207 const Real x = std::exp(mesher_->location(iter, 0));
208 const Real spot = std::min(leverageFct_->maxStrike(),
209 std::max(leverageFct_->minStrike(), x));
210 v[nx] = std::max(0.01, leverageFct_->localVol(time, spot, true));
211 }
212 else {
213 v[iter.index()] = v[nx];
214 }
215 }
216 return v;
217 }
218
219 std::vector<SparseMatrix> FdmHestonFwdOp::toMatrixDecomp() const {
220
221 std::vector<SparseMatrix> retVal(3);
222
223 retVal[0] = mapX_->toMatrix();
224 retVal[1] = mapY_->toMatrix();
225 retVal[2] = correlation_->toMatrix();
226
227 return retVal;
228 }
229
230}
const Real kappa_
const Real correlation_
1-D array used in linear algebra.
Definition: array.hpp:52
const ext::shared_ptr< FirstDerivativeOp > dxMap_
Size size() const override
Array apply_direction(Size direction, const Array &r) const override
const ext::shared_ptr< NinePointLinearOp > correlation_
Array preconditioner(const Array &r, Real s) const override
std::vector< SparseMatrix > toMatrixDecomp() const override
FdmHestonFwdOp(const ext::shared_ptr< FdmMesher > &mesher, const ext::shared_ptr< HestonProcess > &process, FdmSquareRootFwdOp::TransformationType type=FdmSquareRootFwdOp::Plain, ext::shared_ptr< LocalVolTermStructure > leverageFct=ext::shared_ptr< LocalVolTermStructure >(), Real mixingFactor=1.0)
const ext::shared_ptr< LocalVolTermStructure > leverageFct_
void setTime(Time t1, Time t2) override
Time is required.
const ext::shared_ptr< ModTripleBandLinearOp > boundary_
Array apply_mixed(const Array &r) const override
const ext::shared_ptr< YieldTermStructure > qTS_
const ext::shared_ptr< FdmMesher > mesher_
const ext::shared_ptr< TripleBandLinearOp > mapX_
Array getLeverageFctSlice(Time t1, Time t2) const
Array solve_splitting(Size direction, const Array &r, Real s) const override
const ext::shared_ptr< ModTripleBandLinearOp > dxxMap_
const ext::shared_ptr< YieldTermStructure > rTS_
const FdmSquareRootFwdOp::TransformationType type_
Array apply(const Array &r) const override
const ext::shared_ptr< FdmSquareRootFwdOp > mapY_
const DefaultType & t
#define QL_FAIL(message)
throw an error (possibly with file and line information)
Definition: errors.hpp:92
const ext::shared_ptr< YieldTermStructure > rTS_
Heston Fokker-Planck forward operator.
memory layout of a fdm linear operator
mesher for a fdm grid
const ext::shared_ptr< FdmMesher > mesher_
first derivative linear operator
Real Time
continuous quantity with 1-year units
Definition: types.hpp:62
QL_REAL Real
real number
Definition: types.hpp:50
Real Rate
interest rates
Definition: types.hpp:70
std::size_t Size
size of a container
Definition: types.hpp:58
Heston stochastic process.
Real kappa
Real theta
Real v0
Real rho
Real sigma
modifiable triple band linear operator
Definition: any.hpp:35
Array Log(const Array &v)
Definition: array.hpp:861
Array Exp(const Array &v)
Definition: array.hpp:875
STL namespace.
ext::shared_ptr< YieldTermStructure > q
ext::shared_ptr< YieldTermStructure > r
ext::shared_ptr< BlackVolTermStructure > v
Real beta
Definition: sabr.cpp:200
Real alpha
Definition: sabr.cpp:200
second derivative operator
second order mixed derivative linear operator