QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
fdmsquarerootfwdop.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
8 This file is part of QuantLib, a free-software/open-source library
9 for financial quantitative analysts and developers - http://quantlib.org/
10
11 QuantLib is free software: you can redistribute it and/or modify it
12 under the terms of the QuantLib license. You should have received a
13 copy of the license along with this program; if not, please email
14 <quantlib-dev@lists.sf.net>. The license is also available online at
15 <http://quantlib.org/license.shtml>.
16
17 This program is distributed in the hope that it will be useful, but WITHOUT
18 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
19 FOR A PARTICULAR PURPOSE. See the license for more details.
20*/
21
22/*! \file fdmsquarerootfwdop.cpp
23 \brief Fokker-Planck forward operator for an square root process
24*/
25
32
33namespace QuantLib {
34
36 const ext::shared_ptr<FdmMesher>& mesher,
38 Size direction, TransformationType transform)
39 : direction_(direction),
43 transform_(transform),
44 mapX_(transform == Plain ?
46 .mult(kappa*(mesher->locations(direction_)-theta) + sigma*sigma)
47 .add(SecondDerivativeOp(direction_, mesher)
48 .mult(0.5*sigma*sigma*mesher->locations(direction_)))
49 .add(Array(mesher->layout()->size(), kappa)))
50
51 : transform == Power ? new ModTripleBandLinearOp(
53 .mult(0.5*sigma*sigma*mesher->locations(direction_))
54 .add(FirstDerivativeOp(direction_, mesher)
55 .mult(kappa*(mesher->locations(direction_)+theta)))
56 .add(Array(mesher->layout()->size(),
58
60 .mult(Exp(-mesher->locations(direction))
61 *( -0.5*sigma*sigma - kappa*theta) + kappa)
62 .add(SecondDerivativeOp(direction_, mesher)
63 .mult(0.5*sigma*sigma*Exp(-mesher->locations(direction))))
64 .add(kappa*theta*Exp(-mesher->locations(direction))))
65 ),
66 v_ (mesher->layout()->dim()[direction_]) {
67
68 for (const auto& iter : *mesher->layout()) {
69 const Real v = mesher->location(iter, direction_);
70 v_[iter.coordinates()[direction_]] = v;
71 }
72
73 // zero flux boundary condition
74 setLowerBC(mesher);
75 setUpperBC(mesher);
76 }
77
79 const ext::shared_ptr<FdmMesher>& mesher) {
80 const Size n = 1;
81 Real alpha, beta, gamma;
82
83 getCoeff(alpha, beta, gamma, n);
85
86 const Real b = -(h(n-1)+h(n))/zeta(n);
87 const Real c = h(n-1)/zetap(n);
88
89 for (const auto& iter : *mesher->layout()) {
90 if (iter.coordinates()[direction_] == 0) {
91 const Size idx = iter.index();
92 mapX_->diag(idx) = beta + f*b; //*v(n-1);
93 mapX_->upper(idx) = gamma + f*c; //*v(n-1);
94 }
95 }
96 }
97
99 const ext::shared_ptr<FdmMesher>& mesher) {
100 const Size n = v_.size();
101 Real alpha, beta, gamma;
102
103 getCoeff(alpha, beta, gamma, n);
105
106 const Real b = (h(n)+h(n-1))/zeta(n);
107 const Real c = -h(n)/zetam(n);
108
109 for (const auto& iter : *mesher->layout()) {
110 if (iter.coordinates()[direction_] == n-1) {
111 const Size idx = iter.index();
112 mapX_->diag(idx) = beta + f*b; //*v(n+1);
113 mapX_->lower(idx) = alpha + f*c; //*v(n+1);
114 }
115 }
116 }
117
119 if (transform == Plain) {
120 return f0Plain();
121 }
122 else if (transform == Power) {
123 return f0Power();
124 }
125 else if (transform == Log) {
126 return f0Log();
127 }
128 else
129 QL_FAIL("unknown transform");
130 }
131
133 if (transform == Plain) {
134 return f1Plain();
135 }
136 else if (transform == Power) {
137 return f1Power();
138 }
139 else if (transform == Log) {
140 return f1Log();
141 }
142 else
143 QL_FAIL("unknown transform");
144 }
145
147 const Size n = 1;
148 const Real a = -(2*h(n-1)+h(n))/zetam(n);
149 const Real alpha = sigma_*sigma_*v(n)/zetam(n) - mu(n)*h(n)/zetam(n);
150 const Real nu = a*v(n-1) + (2*kappa_*(v(n-1)-theta_) + sigma_*sigma_)
151 /(sigma_*sigma_);
152
153 return alpha/nu*v(n-1);
154 }
155
157 const Size n = v_.size();
158 const Real a = (2*h(n)+h(n-1))/zetap(n);
159 const Real gamma = sigma_*sigma_*v(n)/zetap(n) + mu(n)*h(n-1)/zetap(n);
160 const Real nu = a*v(n+1) + (2*kappa_*(v(n+1)-theta_) + sigma_*sigma_)
161 /(sigma_*sigma_);
162
163 return gamma/nu*v(n+1);
164 }
165
167 const Size n = 1;
168 const Real mu = kappa_*(v(n)+theta_);
169 const Real a = -(2*h(n-1)+h(n))/zetam(n);
170 const Real alpha = sigma_*sigma_*v(n)/zetam(n) - mu*h(n)/zetam(n);
171 const Real nu = a*v(n-1) +2*(kappa_*v(n-1)/(sigma_*sigma_));
172
173 return alpha/nu*v(n-1);
174 }
175
177 const Size n = v_.size();
178 const Real mu = kappa_*(v(n)+theta_);
179 const Real a = (2*h(n)+h(n-1))/zetap(n);
180 const Real gamma = sigma_*sigma_*v(n)/zetap(n) + mu*h(n-1)/zetap(n);
181 const Real nu = a*v(n+1) +2*(kappa_*v(n+1)/(sigma_*sigma_));
182
183 return gamma/nu*v(n+1);
184 }
185
187 const Size n = 1;
188 const Real mu = ((-kappa_*theta_-sigma_*sigma_/2.0)*exp(-v(1))+kappa_);
189 const Real a = -(2*h(n-1)+h(n))/zetam(n);
190 const Real alpha = sigma_*sigma_*exp(-v(n))/zetam(n) - mu*h(n)/zetam(n);
191 const Real nu = a*exp(-v(n-1)) + 2*kappa_*(1-theta_*exp(-v(n-1)))
192 /(sigma_*sigma_);
193
194 return alpha/nu*exp(-v(n-1));
195 }
196
198 const Size n = v_.size();
199 const Real mu = ((-kappa_*theta_-sigma_*sigma_/2.0)*exp(-v(n))+kappa_);
200 const Real a = (2*h(n)+h(n-1))/zetap(n);
201 const Real gamma = sigma_*sigma_*exp(-v(n))/zetap(n) + mu*h(n-1)/zetap(n);
202 const Real nu = a*exp(-v(n+1)) + 2*kappa_*(1-theta_*exp(-v(n+1)))
203 /(sigma_*sigma_);
204
205 return gamma/nu*exp(-v(n+1));
206 }
207
209 if (i > 0 && i <= v_.size()) {
210 return v_[i-1];
211 }
212 else if (i == 0) {
213 if (transform_ == Log) {
214 return 2*v_[0] - v_[1];
215// log(std::max(0.5*exp(v_[0]), exp(v_[0] - 0.01 * (v_[1] - v_[0]))));
216 } else {
217 return std::max(0.5*v_[0], v_[0] - 0.01 * (v_[1] - v_[0]));
218 }
219 }
220 else if (i == v_.size()+1) {
221 return v_.back() + (v_.back() - *(v_.end()-2));
222 }
223 else {
224 QL_FAIL("unknown index");
225 }
226 }
227
229 return v(i+1) - v(i);
230 }
232 return kappa_*(v(i) - theta_) + sigma_*sigma_;
233 }
235 return h(i-1)*(h(i-1)+h(i));
236 }
238 return h(i-1)*h(i);
239 }
241 return h(i)*(h(i-1)+h(i));
242 }
243
245 return 1;
246 }
248 }
249
251 Real& gamma, Size n) const {
252 if (transform_ == Plain) {
253 getCoeffPlain(alpha, beta, gamma, n);
254 }
255 else if (transform_ == Power) {
256 getCoeffPower(alpha, beta, gamma, n);
257 }
258 else if (transform_ == Log) {
259 getCoeffLog(alpha, beta, gamma, n);
260 }
261 }
262
264 Real& gamma, Size n) const {
265 alpha = sigma_*sigma_*v(n)/zetam(n) - mu(n)*h(n)/zetam(n);
266 beta = - sigma_*sigma_*v(n)/zeta(n)
267 + mu(n)*(h(n)-h(n-1))/zeta(n) + kappa_;
268 gamma = sigma_*sigma_*v(n)/zetap(n) + mu(n)*h(n-1)/zetap(n);
269
270 }
271
273 Real& gamma, Size n) const {
274 const Real mu = ((-kappa_*theta_-sigma_*sigma_/2.0)*exp(-v(n))+kappa_);
275 alpha = sigma_*sigma_*exp(-v(n))/zetam(n) - mu*h(n)/zetam(n);
276 beta = - sigma_*sigma_*exp(-v(n))/zeta(n)
277 + mu*(h(n)-h(n-1))/zeta(n) + kappa_*theta_*exp(-v(n));
278 gamma = sigma_*sigma_*exp(-v(n))/zetap(n) + mu*h(n-1)/zetap(n);
279 }
280
282 Real& gamma, Size n) const {
283 const Real mu = kappa_*(theta_+v(n));
284 alpha = (sigma_*sigma_*v(n) - mu*h(n))/zetam(n);
285 beta = (-sigma_*sigma_*v(n) + mu*(h(n)-h(n-1)))/zeta(n)
287 gamma= (sigma_*sigma_*v(n) + mu*h(n-1))/zetap(n);
288 }
289
291 return mapX_->apply(p);
292 }
293
295 return Array(r.size(), 0.0);
296 }
297
299 Size direction, const Array& r) const {
300 if (direction == direction_) {
301 return mapX_->apply(r);
302 }
303 else {
304 return Array(r.size(), 0.0);
305 }
306 }
307
309 Size direction, const Array& r, Real dt) const {
310 if (direction == direction_) {
311 return mapX_->solve_splitting(r, dt, 1.0);
312 }
313 else {
314 return r;
315 }
316 }
317
319 const Array& r, Real dt) const {
320 return solve_splitting(direction_, r, dt);
321 }
322
323 std::vector<SparseMatrix> FdmSquareRootFwdOp::toMatrixDecomp() const {
324 return std::vector<SparseMatrix>(1, mapX_->toMatrix());
325 }
326
327}
const Real kappa_
1-D array used in linear algebra.
Definition: array.hpp:52
const_iterator end() const
Definition: array.hpp:511
Real back() const
Definition: array.hpp:458
Size size() const
dimension of the array
Definition: array.hpp:495
Array apply_direction(Size direction, const Array &r) const override
Array preconditioner(const Array &r, Real s) const override
void setLowerBC(const ext::shared_ptr< FdmMesher > &mesher)
FdmSquareRootFwdOp(const ext::shared_ptr< FdmMesher > &mesher, Real kappa, Real theta, Real sigma, Size direction, TransformationType type=Plain)
std::vector< SparseMatrix > toMatrixDecomp() const override
void getCoeffLog(Real &alpha, Real &beta, Real &gamma, Size n) const
ext::shared_ptr< ModTripleBandLinearOp > mapX_
Real upperBoundaryFactor(TransformationType type=Plain) const
void setUpperBC(const ext::shared_ptr< FdmMesher > &mesher)
void setTime(Time t1, Time t2) override
Time is required.
void getCoeffPlain(Real &alpha, Real &beta, Real &gamma, Size n) const
Array apply_mixed(const Array &r) const override
void getCoeffPower(Real &alpha, Real &beta, Real &gamma, Size n) const
Array solve_splitting(Size direction, const Array &r, Real s) const override
const TransformationType transform_
Array apply(const Array &r) const override
void getCoeff(Real &alpha, Real &beta, Real &gamma, Size n) const
Real lowerBoundaryFactor(TransformationType type=Plain) const
#define QL_FAIL(message)
throw an error (possibly with file and line information)
Definition: errors.hpp:92
ext::function< Real(Real)> b
memory layout of a fdm linear operator
mesher for a fdm grid
Square root linear operator for the Fokker-Planck forward equation.
const Size direction_
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
std::size_t Size
size of a container
Definition: types.hpp:58
Real kappa
Real theta
Real sigma
modifiable triple band linear operator
Definition: any.hpp:35
Array Exp(const Array &v)
Definition: array.hpp:875
ext::shared_ptr< YieldTermStructure > r
ext::shared_ptr< BlackVolTermStructure > v
Real beta
Definition: sabr.cpp:200
Real nu
Definition: sabr.cpp:200
Real alpha
Definition: sabr.cpp:200
second derivative operator