QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
localvolrndcalculator.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) 2015 Johannes Göttker-Schnetmann
5 Copyright (C) 2015 Klaus Spanderen
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 localvolendcalculator.cpp
22 \brief local volatility risk neutral terminal density calculation
23*/
24
35#include <ql/quote.hpp>
38#include <ql/timegrid.hpp>
39#include <utility>
40
41
42namespace QuantLib {
44 ext::shared_ptr<Quote> spot,
45 ext::shared_ptr<YieldTermStructure> rTS,
46 ext::shared_ptr<YieldTermStructure> qTS,
47 const ext::shared_ptr<LocalVolTermStructure>& localVol,
48 Size xGrid,
49 Size tGrid,
50 Real x0Density,
51 Real eps,
52 Size maxIter,
53 Time gaussianStepSize)
54 : xGrid_(xGrid), tGrid_(tGrid), x0Density_(x0Density), localVolProbEps_(eps), maxIter_(maxIter),
55 gaussianStepSize_(gaussianStepSize), spot_(std::move(spot)), localVol_(localVol),
56 rTS_(std::move(rTS)), qTS_(std::move(qTS)),
57 timeGrid_(new TimeGrid(localVol->maxTime(), tGrid)), xm_(tGrid),
58 pm_(new Matrix(tGrid, xGrid)) {
63 }
64
66 ext::shared_ptr<YieldTermStructure> rTS,
67 ext::shared_ptr<YieldTermStructure> qTS,
68 ext::shared_ptr<LocalVolTermStructure> localVol,
69 const ext::shared_ptr<TimeGrid>& timeGrid,
70 Size xGrid,
71 Real x0Density,
72 Real eps,
73 Size maxIter,
74 Time gaussianStepSize)
75 : xGrid_(xGrid), tGrid_(timeGrid->size() - 1), x0Density_(x0Density), localVolProbEps_(eps),
76 maxIter_(maxIter), gaussianStepSize_(gaussianStepSize), spot_(std::move(spot)),
77 localVol_(std::move(localVol)), rTS_(std::move(rTS)), qTS_(std::move(qTS)),
78 timeGrid_(timeGrid), xm_(tGrid_), pm_(new Matrix(tGrid_, xGrid_)) {
83 }
84
86 calculate();
87
88 QL_REQUIRE(t > 0, "positive time expected");
89 QL_REQUIRE(t <= timeGrid_->back(),
90 "given time exceeds local vol time grid");
91
92 const Time tMin = std::min(timeGrid_->at(1), 1.0/365);
93
94 if (t <= tMin) {
95 const Volatility vol = localVol_->localVol(0.0, spot_->value());
96 const Volatility stdDev = vol * std::sqrt(t);
97 const Real xm = - 0.5 * stdDev * stdDev +
98 std::log(spot_->value() * qTS_->discount(t)/rTS_->discount(t));
99
100 return GaussianDistribution(xm, stdDev)(x);
101 }
102 else if (t <= timeGrid_->at(1)) {
103 const Volatility vol = localVol_->localVol(0.0, spot_->value());
104 const Volatility stdDev = vol * std::sqrt(tMin);
105 const Real xm = - 0.5 * stdDev * stdDev +
106 std::log(spot_->value() * qTS_->discount(tMin)/rTS_->discount(tMin));
107
108 const GaussianDistribution gaussianPDF(xm, stdDev);
109
110 const Time deltaT = timeGrid_->at(1) - tMin;
111 return gaussianPDF(x)*(timeGrid_->at(1) - t)/deltaT
112 + probabilityInterpolation(0, x)*(t - tMin)/deltaT;
113 }
114 else {
116 = std::lower_bound(timeGrid_->begin(), timeGrid_->end(), t);
117 const TimeGrid::const_iterator llb = lb-1;
118
119 const Size idx = std::distance(timeGrid_->begin(), lb)-1;
120
121 const Time deltaT = *lb - *llb;
122 return probabilityInterpolation(idx-1, x)*(*lb - t)/deltaT
123 + probabilityInterpolation(idx, x)*(t - *llb)/deltaT;
124 }
125 }
126
128 calculate();
129
130 // get the left side of the integral
131 const Time tc = timeGrid_->closestTime(t);
132 const Size idx = (tc > t) ? timeGrid_->index(tc)-1
133 : std::min(xm_.size()-1, timeGrid_->index(tc));
134
135 Real xl = xm_[idx]->locations().front();
136 Real xr = xm_[idx]->locations().back();
137
138 if (x < xl)
139 return 0.0;
140 else if (x > xr)
141 return 1.0;
142
143 Real addition = 0.1*(xr-xl);
144
145 // left or right hand integral
146 if (x > 0.5*(xr+xl)) {
147 while (pdf(xr, t) > 0.01*localVolProbEps_)
148 {
149 addition*=1.1;
150 xr+=addition;
151 }
152
154 [&](Real _x){ return pdf(_x, t); }, x, xr);
155 }
156 else {
157 while (pdf(xl, t) > 0.01*localVolProbEps_)
158 {
159 addition*=1.1;
160 xl-=addition;
161 }
162
164 [&](Real _x){ return pdf(_x, t); }, xl, x);
165 }
166 }
167
169 calculate();
170
171 const Time closeGridTime(timeGrid_->closestTime(t));
172
173 if (closeGridTime == 0.0) {
174 const Real stepSize = 0.02*(
175 xm_[0]->locations().back() - xm_[0]->locations().front());
177 this, std::log(spot_->value()),
178 0.1*localVolProbEps_, maxIter_, stepSize).inverseCDF(p, t);
179 }
180 else {
181 Array xp(xGrid_);
182 const Size idx = timeGrid_->index(closeGridTime)-1;
183
184 const Array x(xm_[idx]->locations().begin(),
185 xm_[idx]->locations().end());
186 const Real stepSize = 0.005*(x.back() - x.front());
187
188 std::transform(x.begin(), x.end(), pm_->row_begin(idx), xp.begin(), std::multiplies<>());
189
190 const Real xm = DiscreteSimpsonIntegral()(x, xp);
192 this, xm, 0.1*localVolProbEps_, maxIter_, stepSize).inverseCDF(p, t);
193 }
194 }
195
196 ext::shared_ptr<Fdm1dMesher>
198 calculate();
199
200 const Size idx = timeGrid_->index(t);
201 QL_REQUIRE(idx <= xm_.size(), "inconsistent time " << t << " given");
202
203 if (idx > 0) {
204 return xm_[idx-1];
205 }
206 else {
207 return ext::make_shared<Predefined1dMesher>(
208 std::vector<Real>(xGrid_, std::log(spot_->value())));
209 }
210 }
211
212 ext::shared_ptr<TimeGrid> LocalVolRNDCalculator::timeGrid() const {
213 return timeGrid_;
214 }
215
217 rescaleTimeSteps_.clear();
218
219 const Time sT = timeGrid_->at(1);
220 Time t = std::min(sT, (gaussianStepSize_ > 0.0) ? gaussianStepSize_
221 : 0.5*sT);
222 const Volatility vol = localVol_->localVol(0.0, spot_->value());
223
224 const Volatility stdDev = vol * std::sqrt(t);
225 Real xm = - 0.5 * stdDev * stdDev +
226 std::log(spot_->value() * qTS_->discount(t)/rTS_->discount(t));
227
228 const Volatility stdDevOfFirstStep = vol * std::sqrt(sT);
229 const Real normInvEps = InverseCumulativeNormal()(1 - localVolProbEps_);
230
231 Real sLowerBound = xm - normInvEps * stdDevOfFirstStep;
232 Real sUpperBound = xm + normInvEps * stdDevOfFirstStep;
233
234 ext::shared_ptr<Fdm1dMesher> mesher(
235 new Concentrating1dMesher(sLowerBound, sUpperBound, xGrid_,
236 std::make_pair(xm, x0Density_), true));
237
238 Array p(mesher->size());
239 Array x(mesher->locations().begin(), mesher->locations().end());
240
241 const GaussianDistribution gaussianPDF(xm, vol * std::sqrt(t));
242
243 for (Size idx=0; idx < p.size(); ++idx) {
244 p[idx] = gaussianPDF(x[idx]);
245 }
246 p = rescalePDF(x, p);
247
248 QL_REQUIRE(x.size() > 10, "x grid is too small. "
249 "Minimum size is greater than 10");
250
251 const Size b = std::max(Size(1), Size(x.size()*0.04));
252
253 ext::shared_ptr<DouglasScheme> evolver(
254 new DouglasScheme(0.5,
255 ext::make_shared<FdmLocalVolFwdOp>(
256 ext::make_shared<FdmMesherComposite>(mesher),
257 spot_, rTS_, qTS_, localVol_)));
258
259 pFct_.resize(tGrid_);
260
261 for (Size i=1; i <= tGrid_; ++i) {
262 const Time dt = timeGrid_->at(i) - t;
263
264 // leaking probability mass?
265 const Real maxLeftValue =
266 std::max(std::fabs(*std::min_element(p.begin(), p.begin()+b)),
267 std::fabs(*std::max_element(p.begin(), p.begin()+b)));
268 const Real maxRightValue =
269 std::max(std::fabs(*std::min_element(p.end()-b, p.end())),
270 std::fabs(*std::max_element(p.end()-b, p.end())));
271
272 if (std::max(maxLeftValue, maxRightValue) > localVolProbEps_) {
273 rescaleTimeSteps_.push_back(i);
274
275 const Real oldLowerBound = sLowerBound;
276 const Real oldUpperBound = sUpperBound;
277
278 xm = DiscreteSimpsonIntegral()(x, x*p);
279 Array vols(x.size());
280 for (Size j=0; j < vols.size(); ++j) {
281 vols[j] = localVol_->localVol(t + dt, std::exp(x[j]));
282 }
283
284 const Real vm = DiscreteSimpsonIntegral()(x, vols)
285 /(x.back() - x.front());
286
287 const Real scalingFactor = vm*std::sqrt(0.5*timeGrid_->back());
288
289 if (maxLeftValue > localVolProbEps_)
290 sLowerBound -= scalingFactor*(oldUpperBound-oldLowerBound);
291 if (maxRightValue > localVolProbEps_)
292 sUpperBound += scalingFactor*(oldUpperBound-oldLowerBound);
293
294 mesher = ext::shared_ptr<Fdm1dMesher>(
295 new Concentrating1dMesher(sLowerBound, sUpperBound, xGrid_,
296 std::make_pair(xm, 0.1), false));
297
298 const CubicNaturalSpline pSpline(x.begin(), x.end(), p.begin());
299 const Array xn(mesher->locations().begin(),
300 mesher->locations().end());
301 Array pn(xn.size(), 0.0);
302
303 for (Size j=0; j < xn.size(); ++j) {
304 if (xn[j] >= oldLowerBound && xn[j] <= oldUpperBound)
305 pn[j] = pSpline(xn[j]);
306 }
307
308 x = xn;
309 p = rescalePDF(xn, pn);
310
311 evolver = ext::make_shared<DouglasScheme>(0.5,
312 ext::make_shared<FdmLocalVolFwdOp>(
313 ext::make_shared<FdmMesherComposite>(mesher),
315 }
316 evolver->setStep(dt);
317 t+=dt;
318
319 if (dt > QL_EPSILON) {
320 evolver->step(p, t);
321 p = rescalePDF(x, p);
322 }
323
324 xm_[i-1] = mesher;
325 std::copy(p.begin(), p.end(), pm_->row_begin(i-1));
326 pFct_[i-1] = ext::make_shared<CubicNaturalSpline>(
327 xm_[i-1]->locations().begin(),
328 xm_[i-1]->locations().end(),
329 pm_->row_begin(i-1));
330 }
331 }
332
333
334 std::vector<Size> LocalVolRNDCalculator::rescaleTimeSteps() const {
335 calculate();
336
337 return rescaleTimeSteps_;
338 }
339
341 Size idx, Real x) const {
342 calculate();
343
344 if ( x < xm_[idx]->locations().front()
345 || x > xm_[idx]->locations().back())
346 return 0.0;
347 else
348 return (*pFct_[idx])(x);
349 }
350
352 return p/DiscreteSimpsonIntegral()(x, p);
353 }
354
355}
356
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
Real front() const
Definition: array.hpp:451
const_iterator begin() const
Definition: array.hpp:503
Integral of a one-dimensional function.
Inverse cumulative normal distribution function.
virtual void calculate() const
Definition: lazyobject.hpp:253
std::vector< ext::shared_ptr< Interpolation > > pFct_
const ext::shared_ptr< Quote > spot_
LocalVolRNDCalculator(ext::shared_ptr< Quote > spot, ext::shared_ptr< YieldTermStructure > rTS, ext::shared_ptr< YieldTermStructure > qTS, const ext::shared_ptr< LocalVolTermStructure > &localVol, Size xGrid=101, Size tGrid=51, Real x0Density=0.1, Real localVolProbEps=1e-6, Size maxIter=10000, Time gaussianStepSize=-Null< Time >())
Array rescalePDF(const Array &x, const Array &p) const
const ext::shared_ptr< YieldTermStructure > qTS_
std::vector< Size > rescaleTimeSteps() const
ext::shared_ptr< TimeGrid > timeGrid() const
Real pdf(Real x, Time t) const override
Real cdf(Real x, Time t) const override
const ext::shared_ptr< YieldTermStructure > rTS_
const ext::shared_ptr< TimeGrid > timeGrid_
const ext::shared_ptr< Matrix > pm_
ext::shared_ptr< Fdm1dMesher > mesher(Time t) const
std::vector< ext::shared_ptr< Fdm1dMesher > > xm_
Real invcdf(Real p, Time t) const override
Real probabilityInterpolation(Size idx, Real x) const
const ext::shared_ptr< LocalVolTermStructure > localVol_
Matrix used in linear algebra.
Definition: matrix.hpp:41
Normal distribution function.
std::pair< iterator, bool > registerWith(const ext::shared_ptr< Observable > &)
Definition: observable.hpp:228
time grid class
Definition: timegrid.hpp:43
std::vector< Time >::const_iterator const_iterator
Definition: timegrid.hpp:158
One-dimensional grid mesher concentrating around critical points.
cubic interpolation between discrete points
const DefaultType & t
integrals on non uniform grids
Douglas operator splitting.
#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
const ext::shared_ptr< YieldTermStructure > rTS_
local volatility linear operator for the Fokker-Planck forward equation
FdmMesher which is a composite of Fdm1dMesher.
integral of a one-dimensional function using the adaptive Gauss-Lobatto integral
#define QL_EPSILON
Definition: qldefines.hpp:178
Real Time
continuous quantity with 1-year units
Definition: types.hpp:62
QL_REAL Real
real number
Definition: types.hpp:50
Real Volatility
volatility
Definition: types.hpp:78
std::size_t Size
size of a container
Definition: types.hpp:58
local volatility risk neutral terminal density calculation
Local volatility term structure base class.
Definition: any.hpp:35
NormalDistribution GaussianDistribution
STL namespace.
normal, cumulative and inverse cumulative distributions
One-dimensional mesher build from a given set of points.
purely virtual base class for market observables
discrete time grid
Interest-rate term structure.