1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
|
/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
/*
Copyright (C) 2014 Klaus Spanderen
This file is part of QuantLib, a free-software/open-source library
for financial quantitative analysts and developers - http://quantlib.org/
QuantLib is free software: you can redistribute it and/or modify it
under the terms of the QuantLib license. You should have received a
copy of the license along with this program; if not, please email
<quantlib-dev@lists.sf.net>. The license is also available online at
<https://www.quantlib.org/license.shtml>.
This program is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
FOR A PARTICULAR PURPOSE. See the license for more details.
*/
/*! \file fdmhestongreenfct.cpp
\brief Heston Fokker-Planck Green's function
*/
#include <ql/math/functional.hpp>
#include <ql/methods/finitedifferences/meshers/fdmmesher.hpp>
#include <ql/methods/finitedifferences/operators/fdmlinearopiterator.hpp>
#include <ql/methods/finitedifferences/operators/fdmlinearoplayout.hpp>
#include <ql/methods/finitedifferences/utilities/fdmhestongreensfct.hpp>
#include <ql/methods/finitedifferences/utilities/squarerootprocessrndcalculator.hpp>
#include <ql/processes/hestonprocess.hpp>
#include <utility>
namespace QuantLib {
FdmHestonGreensFct::FdmHestonGreensFct(ext::shared_ptr<FdmMesher> mesher,
ext::shared_ptr<HestonProcess> process,
FdmSquareRootFwdOp::TransformationType trafoType_,
const Real l0)
: l0_(l0), mesher_(std::move(mesher)), process_(std::move(process)), trafoType_(trafoType_) {}
Array FdmHestonGreensFct::get(Time t, Algorithm algorithm) const {
const Rate r = process_->riskFreeRate()->forwardRate(0, t, Continuous);
const Rate q = process_->dividendYield()->forwardRate(0,t, Continuous);
const Real s0 = process_->s0()->value();
const Real v0 = process_->v0();
const Real x0 = std::log(s0) + (r-q-0.5*v0*l0_*l0_)*t;
const Real rho = process_->rho();
const Real theta = process_->theta();
const Real kappa = process_->kappa();
const Real sigma = process_->sigma();
Array p(mesher_->layout()->size());
for (const auto& iter : *mesher_->layout()) {
const Real x = mesher_->location(iter, 0);
const Real v = (trafoType_ != FdmSquareRootFwdOp::Log)
? mesher_->location(iter, 1)
: std::exp(mesher_->location(iter, 1));
Real retVal;
switch (algorithm) {
case ZeroCorrelation:
{
const Real sd_x = l0_*std::sqrt(v0*t);
const Real p_x = M_1_SQRTPI*M_SQRT1_2/sd_x
* std::exp(-0.5*squared((x - x0)/sd_x));
const Real p_v = SquareRootProcessRNDCalculator(
v0, kappa, theta, sigma).pdf(v, t);
retVal = p_v*p_x;
}
break;
case SemiAnalytical:
retVal = process_->pdf(x, v, t, 1e-4);
break;
case Gaussian:
{
const Real sd_x = l0_*std::sqrt(v0*t);
const Real sd_v = sigma*std::sqrt(v0*t);
const Real z0 = v0 + kappa*(theta - v0)*t;
retVal = 1.0/(M_TWOPI*sd_x*sd_v*std::sqrt(1-rho*rho))
*std::exp(-( squared((x-x0)/sd_x)
+ squared((v-z0)/sd_v)
- 2*rho*(x-x0)*(v-z0)/(sd_x*sd_v))
/(2*(1-rho*rho)) );
}
break;
default:
QL_FAIL("unknown algorithm");
}
switch (trafoType_) {
case FdmSquareRootFwdOp::Log:
retVal*=v;
break;
case FdmSquareRootFwdOp::Plain:
break;
case FdmSquareRootFwdOp::Power:
retVal*=std::pow(v, 1.0 - 2*kappa*theta/(sigma*sigma));
break;
default:
QL_FAIL("unknown transformation type");
}
p[iter.index()] = retVal;
}
return p;
}
}
|
