QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
fdmhestonvariancemesher.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) 2008 Andreas Gaida
5 Copyright (C) 2008 Ralph Schreyer
6 Copyright (C) 2008, 2019 Klaus Spanderen
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
29#include <boost/accumulators/accumulators.hpp>
30#include <boost/accumulators/statistics/mean.hpp>
31#include <boost/accumulators/statistics/stats.hpp>
32#include <set>
33#include <algorithm>
34
35namespace QuantLib {
36
37 namespace {
38 struct interpolated_volatility {
39 interpolated_volatility(const std::vector<Real>& pGrid,
40 const std::vector<Real>& vGrid)
41 : variance(pGrid.begin(), pGrid.end(), vGrid.begin()) {}
42 Real operator()(Real x) const {
43 return std::sqrt(variance(x, true));
44 }
45 LinearInterpolation variance;
46 };
47 }
48
50 Size size,
51 const ext::shared_ptr<HestonProcess> & process,
52 Time maturity, Size tAvgSteps, Real epsilon,
53 Real mixingFactor)
54 : Fdm1dMesher(size) {
55
56 std::vector<Real> vGrid(size, 0.0), pGrid(size, 0.0);
57 const Real mixedSigma = process->sigma()*mixingFactor;
58 const Real df = 4*process->theta()*process->kappa()/squared(mixedSigma);
59 try {
60 std::multiset<std::pair<Real, Real> > grid;
61
62 for (Size l=1; l<=tAvgSteps; ++l) {
63 const Real t = (maturity*l)/tAvgSteps;
64 const Real ncp = 4*process->kappa()*std::exp(-process->kappa()*t)/(squared(mixedSigma)
65 *(1-std::exp(-process->kappa()*t)))*process->v0();
66 const Real k = squared(mixedSigma)
67 *(1-std::exp(-process->kappa()*t))/(4*process->kappa());
68
69 const Real qMin = 0.0; // v_min = 0.0;
70 const Real qMax = std::max(process->v0(),
72 df, ncp, 100, 1e-8)(1-epsilon));
73
74 const Real minVStep=(qMax-qMin)/(50*size);
75 Real ps,p = 0.0;
76
77 Real vTmp = qMin;
78 grid.insert(std::pair<Real, Real>(qMin, epsilon));
79
80 for (Size i=1; i < size; ++i) {
81 ps = (1 - epsilon - p)/(size-i);
82 p += ps;
84 df, ncp, 100, 1e-8)(p);
85
86 const Real vx = std::max(vTmp+minVStep, tmp);
88 vTmp=vx;
89 grid.insert(std::pair<Real, Real>(vx, p));
90 }
91 }
92 QL_REQUIRE(grid.size() == size*tAvgSteps,
93 "something wrong with the grid size");
94
95 const std::vector<std::pair<Real, Real> > tp(grid.begin(), grid.end());
96
97 for (Size i=0; i < size; ++i) {
98 const Size b = (i*tp.size())/size;
99 const Size e = ((i+1)*tp.size())/size;
100 for (Size j=b; j < e; ++j) {
101 vGrid[i]+=tp[j].first/(e-b);
102 pGrid[i]+=tp[j].second/(e-b);
103 }
104 }
105 }
106 catch (const Error&) {
107 // use default mesh
108 const Real vol = mixedSigma*
109 std::sqrt(process->theta()/(2*process->kappa()));
110
111 const Real mean = process->theta();
112 const Real upperBound = std::max(process->v0()+4*vol, mean+4*vol);
113 const Real lowerBound
114 = std::max(0.0, std::min(process->v0()-4*vol, mean-4*vol));
115
116 for (Size i=0; i < size; ++i) {
117 pGrid[i] = i/(size-1.0);
118 vGrid[i] = lowerBound + i*(upperBound-lowerBound)/(size-1.0);
119 }
120 }
121
122 Real skewHint = ((process->kappa() != 0.0)
123 ? Real(std::max(1.0, mixedSigma/process->kappa())) : 1.0);
124
125 std::sort(pGrid.begin(), pGrid.end());
126 volaEstimate_ = GaussLobattoIntegral(100000, 1e-4)(
127 interpolated_volatility(pGrid, vGrid),
128 pGrid.front(), pGrid.back())*std::pow(skewHint, 1.5);
129
130 const Real v0 = process->v0();
131 for (Size i=1; i<vGrid.size(); ++i) {
132 if (vGrid[i-1] <= v0 && vGrid[i] >= v0) {
133 if (std::fabs(vGrid[i-1] - v0) < std::fabs(vGrid[i] - v0))
134 vGrid[i-1] = v0;
135 else
136 vGrid[i] = v0;
137 }
138 }
139
140 std::copy(vGrid.begin(), vGrid.end(), locations_.begin());
141
142 for (Size i=0; i < size-1; ++i) {
143 dminus_[i+1] = dplus_[i] = vGrid[i+1] - vGrid[i];
144 }
145 dplus_.back() = dminus_.front() = Null<Real>();
146 }
147
148
150 Size size,
151 const ext::shared_ptr<HestonProcess>& process,
152 const ext::shared_ptr<LocalVolTermStructure>& leverageFct,
153 Time maturity, Size tAvgSteps, Real epsilon,
154 Real mixingFactor)
155 : Fdm1dMesher(size) {
156
157 const FdmHestonVarianceMesher mesher(
158 size, process, maturity, tAvgSteps, epsilon, mixingFactor);
159
160 for (Size i=0; i < size; ++i) {
161 dplus_[i] = mesher.dplus(i);
162 dminus_[i] = mesher.dminus(i);
163 locations_[i] = mesher.location(i);
164 }
165
166 volaEstimate_ = mesher.volaEstimate();
167
168 if (leverageFct != nullptr) {
169 typedef boost::accumulators::accumulator_set<
170 Real, boost::accumulators::stats<
171 boost::accumulators::tag::mean> >
172 accumulator_set;
173
174 accumulator_set acc;
175
176 const Real s0 = process->s0()->value();
177
178 acc(leverageFct->localVol(0, s0, true));
179
180 const Handle<YieldTermStructure> rTS = process->riskFreeRate();
181 const Handle<YieldTermStructure> qTS = process->dividendYield();
182
183 for (Size l=1; l <= tAvgSteps; ++l) {
184 const Real t = (maturity*l)/tAvgSteps;
185 const Real vol = volaEstimate_ * boost::accumulators::mean(acc);
186
187 const Real fwd = s0*qTS->discount(t)/rTS->discount(t);
188
189 const Size sAvgSteps = 50;
190
191 std::vector<Real> u(sAvgSteps), sig(sAvgSteps);
192
193 for (Size i=0; i < sAvgSteps; ++i) {
194 u[i] = epsilon + ((1-2*epsilon)/(sAvgSteps-1))*i;
195 const Real x = InverseCumulativeNormal()(u[i]);
196
197 const Real gf = x*vol*std::sqrt(t);
198 const Real f = fwd*std::exp(gf);
199
200 sig[i] = squared(leverageFct->localVol(t, f, true));
201 }
202
203 const Real leverageAvg =
204 GaussLobattoIntegral(10000, 1e-4)(
205 interpolated_volatility(u, sig), u.front(), u.back())
206 / (1-2*epsilon);
207
208 acc(leverageAvg);
209 }
210 volaEstimate_ *= boost::accumulators::mean(acc);
211 }
212 }
213}
Chi-square (central and non-central) distributions.
Base error class.
Definition: errors.hpp:39
std::vector< Real > locations_
Definition: fdm1dmesher.hpp:47
Real dminus(Size index) const
Definition: fdm1dmesher.hpp:42
Real location(Size index) const
Definition: fdm1dmesher.hpp:43
std::vector< Real > dplus_
Definition: fdm1dmesher.hpp:48
std::vector< Real > dminus_
Definition: fdm1dmesher.hpp:48
Real dplus(Size index) const
Definition: fdm1dmesher.hpp:41
FdmHestonLocalVolatilityVarianceMesher(Size size, const ext::shared_ptr< HestonProcess > &process, const ext::shared_ptr< LocalVolTermStructure > &leverageFct, Time maturity, Size tAvgSteps=10, Real epsilon=0.0001, Real mixingFactor=1.0)
FdmHestonVarianceMesher(Size size, const ext::shared_ptr< HestonProcess > &process, Time maturity, Size tAvgSteps=10, Real epsilon=0.0001, Real mixingFactor=1.0)
Integral of a one-dimensional function.
Shared handle to an observable.
Definition: handle.hpp:41
Inverse cumulative normal distribution function.
template class providing a null value for a given type.
Definition: null.hpp:76
const DefaultType & t
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
Definition: errors.hpp:117
ext::function< Real(Real)> b
LinearInterpolation variance
One-dimensional grid mesher for the variance part of the Heston model.
integral of a one-dimensional function using the adaptive Gauss-Lobatto integral
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 v0
linear interpolation between discrete points
Local volatility term structure base class.
functionals and combinators not included in the STL
Definition: any.hpp:35
T squared(T x)
Definition: functional.hpp:37
normal, cumulative and inverse cumulative distributions