22#include <ql/models/equity/hestonslvfdmmodel.hpp>
23#include <ql/math/distributions/normaldistribution.hpp>
24#include <ql/math/integrals/discreteintegrals.hpp>
25#include <ql/math/interpolations/bilinearinterpolation.hpp>
26#include <ql/methods/finitedifferences/meshers/concentrating1dmesher.hpp>
27#include <ql/methods/finitedifferences/meshers/fdmmeshercomposite.hpp>
28#include <ql/methods/finitedifferences/meshers/predefined1dmesher.hpp>
29#include <ql/methods/finitedifferences/operators/fdmlinearoplayout.hpp>
30#include <ql/methods/finitedifferences/operators/fdmhestonfwdop.hpp>
31#include <ql/methods/finitedifferences/schemes/craigsneydscheme.hpp>
32#include <ql/methods/finitedifferences/schemes/douglasscheme.hpp>
33#include <ql/methods/finitedifferences/schemes/expliciteulerscheme.hpp>
34#include <ql/methods/finitedifferences/schemes/hundsdorferscheme.hpp>
35#include <ql/methods/finitedifferences/schemes/impliciteulerscheme.hpp>
36#include <ql/methods/finitedifferences/schemes/modifiedcraigsneydscheme.hpp>
37#include <ql/methods/finitedifferences/solvers/fdmbackwardsolver.hpp>
38#include <ql/methods/finitedifferences/utilities/fdmmesherintegral.hpp>
39#include <ql/methods/finitedifferences/utilities/localvolrndcalculator.hpp>
40#include <ql/methods/finitedifferences/utilities/squarerootprocessrndcalculator.hpp>
41#include <ql/models/equity/hestonmodel.hpp>
42#include <ql/quotes/simplequote.hpp>
43#include <ql/termstructures/volatility/equityfx/fixedlocalvolsurface.hpp>
44#include <ql/termstructures/volatility/equityfx/localvoltermstructure.hpp>
45#include <ql/timegrid.hpp>
53 ext::shared_ptr<Fdm1dMesher> varianceMesher(
54 const SquareRootProcessRNDCalculator& rnd,
56 Real v0,
const HestonSLVFokkerPlanckFdmParams& params) {
58 std::vector<ext::tuple<Real, Real, bool> > cPoints;
60 const Real v0Density = params.v0Density;
61 const Real upperBoundDensity = params.vUpperBoundDensity;
62 const Real lowerBoundDensity = params.vLowerBoundDensity;
64 Real lowerBound = Null<Real>(), upperBound = -Null<Real>();
66 for (
Size i=0; i <= 10; ++i) {
67 const Time t = t0 + i/10.0*(t1-t0);
68 lowerBound = std::min(
69 lowerBound, rnd.invcdf(params.vLowerEps, t));
70 upperBound = std::max(
71 upperBound, rnd.invcdf(1.0-params.vUpperEps, t));
74 lowerBound = std::max(lowerBound, params.vMin);
75 switch (params.trafoType) {
78 lowerBound = std::log(lowerBound);
79 upperBound = std::log(upperBound);
81 const Real v0Center = std::log(v0);
84 ext::make_tuple(lowerBound, lowerBoundDensity,
false),
85 ext::make_tuple(v0Center, v0Density,
true),
86 ext::make_tuple(upperBound, upperBoundDensity,
false)
89 return ext::make_shared<Concentrating1dMesher>(
90 lowerBound, upperBound, vGrid, cPoints, 1e-8);
95 const Real v0Center = v0;
98 ext::make_tuple(lowerBound, lowerBoundDensity,
false),
99 ext::make_tuple(v0Center, v0Density,
true),
100 ext::make_tuple(upperBound, upperBoundDensity,
false)
103 return ext::make_shared<Concentrating1dMesher>(
104 lowerBound, upperBound, vGrid, cPoints, 1e-8);
109 const Real v0Center = v0;
112 ext::make_tuple(lowerBound, lowerBoundDensity,
false),
113 ext::make_tuple(v0Center, v0Density,
true),
114 ext::make_tuple(upperBound, upperBoundDensity,
false)
117 return ext::make_shared<Concentrating1dMesher>(
118 lowerBound, upperBound, vGrid, cPoints, 1e-8);
122 QL_FAIL(
"transformation type is not implemented");
126 Real integratePDF(
const Array& p,
127 const ext::shared_ptr<FdmMesherComposite>& mesher,
132 return FdmMesherIntegral(
133 mesher, DiscreteSimpsonIntegral()).integrate(p);
137 for (
const auto& iter : *mesher->layout()) {
138 const Size idx = iter.index();
139 const Real nu = mesher->location(iter, 1);
141 tp[idx] = p[idx]*std::pow(nu, alpha-1);
144 return FdmMesherIntegral(
145 mesher, DiscreteSimpsonIntegral()).integrate(tp);
152 const ext::shared_ptr<FdmMesherComposite>& mesher,
155 return p/integratePDF(p, mesher, trafoType, alpha);
159 template <
class Interpolator>
162 const ext::shared_ptr<FdmMesherComposite>& oldMesher,
163 const ext::shared_ptr<FdmMesherComposite>& newMesher,
164 const Interpolator& interp = Interpolator()) {
166 QL_REQUIRE( oldMesher->layout()->size() == newMesher->layout()->size()
167 && oldMesher->layout()->size() == p.size(),
168 "inconsistent mesher or vector size given");
170 Matrix m(oldMesher->layout()->dim()[1], oldMesher->layout()->dim()[0]);
171 for (
Size i=0; i < m.rows(); ++i) {
172 std::copy(p.begin() + i*m.columns(),
173 p.begin() + (i+1)*m.columns(), m.row_begin(i));
175 const Interpolation2D interpol = interp.interpolate(
176 oldMesher->getFdm1dMeshers()[0]->locations().begin(),
177 oldMesher->getFdm1dMeshers()[0]->locations().end(),
178 oldMesher->getFdm1dMeshers()[1]->locations().begin(),
179 oldMesher->getFdm1dMeshers()[1]->locations().end(), m);
181 Array pNew(p.size());
182 for (
const auto& iter : *newMesher->layout()) {
183 const Real x = newMesher->location(iter, 0);
184 const Real v = newMesher->location(iter, 1);
186 if ( x > interpol.xMax() || x < interpol.xMin()
187 || v > interpol.yMax() || v < interpol.yMin() ) {
188 pNew[iter.index()] = 0;
191 pNew[iter.index()] = interpol(x, v);
200 virtual ~FdmScheme() =
default;
201 virtual void step(Array& a,
Time t) = 0;
202 virtual void setStep(
Time dt) = 0;
206 class FdmSchemeWrapper :
public FdmScheme {
208 explicit FdmSchemeWrapper(T* scheme)
209 : scheme_(scheme) { }
211 void step(Array& a,
Time t)
override { scheme_->step(a, t); }
212 void setStep(
Time dt)
override { scheme_->setStep(dt); }
215 const std::unique_ptr<T> scheme_;
218 ext::shared_ptr<FdmScheme> fdmSchemeFactory(
219 const FdmSchemeDesc desc,
220 const ext::shared_ptr<FdmLinearOpComposite>& op) {
224 return ext::shared_ptr<FdmScheme>(
225 new FdmSchemeWrapper<HundsdorferScheme>(
226 new HundsdorferScheme(desc.theta, desc.mu, op)));
228 return ext::shared_ptr<FdmScheme>(
229 new FdmSchemeWrapper<DouglasScheme>(
230 new DouglasScheme(desc.theta, op)));
232 return ext::shared_ptr<FdmScheme>(
233 new FdmSchemeWrapper<CraigSneydScheme>(
234 new CraigSneydScheme(desc.theta, desc.mu, op)));
236 return ext::shared_ptr<FdmScheme>(
237 new FdmSchemeWrapper<ModifiedCraigSneydScheme>(
238 new ModifiedCraigSneydScheme(
239 desc.theta, desc.mu, op)));
241 return ext::shared_ptr<FdmScheme>(
242 new FdmSchemeWrapper<ImplicitEulerScheme>(
243 new ImplicitEulerScheme(op)));
245 return ext::shared_ptr<FdmScheme>(
246 new FdmSchemeWrapper<ExplicitEulerScheme>(
247 new ExplicitEulerScheme(op)));
249 QL_FAIL(
"Unknown scheme type");
259 std::vector<Date> mandatoryDates,
260 const Real mixingFactor)
261 : localVol_(
std::move(localVol)), hestonModel_(
std::move(hestonModel)), endDate_(endDate),
262 params_(
std::move(params)), mandatoryDates_(
std::move(mandatoryDates)),
263 mixingFactor_(mixingFactor), logging_(logging) {
277 ext::shared_ptr<LocalVolTermStructure>
289 const ext::shared_ptr<Quote> spot
291 const ext::shared_ptr<YieldTermStructure> rTS
293 const ext::shared_ptr<YieldTermStructure> qTS
301 const Real alpha = 2*kappa*theta/(mixedSigma*mixedSigma);
307 const Date referenceDate = rTS->referenceDate();
311 QL_REQUIRE(referenceDate <
endDate_,
312 "reference date must be smaller than final calibration date");
315 "final calibration maturity exceeds local volatility surface");
322 std::vector<Time> times(1, tIdx);
326 const Time dt = maxDt*decayFactor + minDt*(1.0-decayFactor);
328 times.push_back(std::min(T, tIdx+=dt));
332 times.push_back(dc.
yearFraction(referenceDate, mandatoryDate));
335 const ext::shared_ptr<TimeGrid> timeGrid(
336 new TimeGrid(times.begin(), times.end()));
346 const std::vector<Size> rescaleSteps
350 v0, kappa, theta, mixedSigma);
355 std::vector<ext::shared_ptr<Fdm1dMesher> > xMesher, vMesher;
356 xMesher.reserve(timeGrid->size());
357 vMesher.reserve(timeGrid->size());
359 xMesher.push_back(localVolRND.
mesher(0.0));
360 vMesher.push_back(ext::make_shared<Predefined1dMesher>(
361 std::vector<Real>(vGrid, v0)));
364 for (
Size i=1; i < timeGrid->size(); ++i) {
365 xMesher.push_back(localVolRND.
mesher(timeGrid->at(i)));
367 if ((rescaleIdx < rescaleSteps.size())
368 && (i == rescaleSteps[rescaleIdx])) {
370 vMesher.push_back(varianceMesher(squareRootRnd,
371 timeGrid->at(rescaleSteps[rescaleIdx-1]),
372 (rescaleIdx < rescaleSteps.size())
373 ? timeGrid->at(rescaleSteps[rescaleIdx])
378 vMesher.push_back(vMesher.back());
382 ext::shared_ptr<FdmMesherComposite> mesher
383 = ext::make_shared<FdmMesherComposite>(
384 xMesher.at(1), vMesher.at(1));
387 =
localVol_->localVol(0.0, spot->value())/std::sqrt(v0);
389 ext::shared_ptr<Matrix> L(
new Matrix(xGrid, timeGrid->size()));
392 std::fill(L->column_begin(0),L->column_end(0), l0);
393 std::fill(L->column_begin(1),L->column_end(1), l0);
396 std::vector<ext::shared_ptr<std::vector<Real> > > vStrikes(
399 for (
Size i=0; i < timeGrid->size(); ++i) {
400 vStrikes[i] = ext::make_shared<std::vector<Real> >(xGrid);
401 if (xMesher[i]->locations().front()
402 == xMesher[i]->locations().back()) {
403 std::fill(vStrikes[i]->begin(), vStrikes[i]->end(),
404 std::exp(xMesher[i]->locations().front()));
407 std::transform(xMesher[i]->locations().begin(),
408 xMesher[i]->locations().end(),
409 vStrikes[i]->begin(),
410 [](
Real x) ->
Real {
return std::exp(x); });
414 const ext::shared_ptr<FixedLocalVolSurface> leverageFct(
417 ext::shared_ptr<FdmLinearOpComposite> hestonFwdOp(
424 const LogEntry entry = { timeGrid->at(1),
425 ext::make_shared<Array>(p), mesher };
429 for (
Size i=2; i < times.size(); ++i) {
430 const Time t = timeGrid->at(i);
431 const Time dt = t - timeGrid->at(i-1);
433 if ( mesher->getFdm1dMeshers()[0] != xMesher[i]
434 || mesher->getFdm1dMeshers()[1] != vMesher[i]) {
435 const ext::shared_ptr<FdmMesherComposite> newMesher(
438 p = reshapePDF<Bilinear>(p, mesher, newMesher);
441 p = rescalePDF(p, mesher, trafoType, alpha);
443 hestonFwdOp = ext::shared_ptr<FdmLinearOpComposite>(
450 Array(mesher->getFdm1dMeshers()[0]->locations().begin(),
451 mesher->getFdm1dMeshers()[0]->locations().end())));
453 mesher->getFdm1dMeshers()[1]->locations().begin(),
454 mesher->getFdm1dMeshers()[1]->locations().end());
463 const ext::shared_ptr<FdmScheme> fdmScheme(
464 fdmSchemeFactory(fdmSchemeDesc, hestonFwdOp));
468 for (
Size k=0; k < vGrid; ++k)
469 pSlice[k] = pn[j + k*xGrid];
481 const Real scale = pInt/vpInt;
484 const Real l = (scale >= 0.0)
487 (*L)[j][i] = std::min(50.0, std::max(0.001, l));
489 leverageFct->setInterpolation(
Linear());
492 const Real sLowerBound = std::max(x.
front(),
493 std::exp(localVolRND.
invcdf(
495 const Real sUpperBound = std::min(x.
back(),
496 std::exp(localVolRND.
invcdf(
499 const Real lowerL = leverageFct->localVol(t, sLowerBound);
500 const Real upperL = leverageFct->localVol(t, sUpperBound);
503 if (x[j] < sLowerBound)
504 std::fill(L->row_begin(j)+i,
505 std::min(L->row_begin(j)+i+1, L->row_end(j)),
507 else if (x[j] > sUpperBound)
508 std::fill(L->row_begin(j)+i,
509 std::min(L->row_begin(j)+i+1, L->row_end(j)),
512 QL_FAIL(
"internal error");
514 leverageFct->setInterpolation(
Linear());
518 fdmScheme->setStep(dt);
519 fdmScheme->step(pn, t);
522 p = rescalePDF(p, mesher, trafoType, alpha);
526 = { t, ext::make_shared<Array>(p), mesher };
1-D array used in linear algebra.
Size size() const
dimension of the array
Time yearFraction(const Date &, const Date &, const Date &refPeriodStart=Date(), const Date &refPeriodEnd=Date()) const
Returns the period between two dates as a fraction of year.
Array get(Time t, Algorithm algorithm) const
Shared handle to an observable.
ext::shared_ptr< LocalVolTermStructure > leverageFunction() const
void performCalculations() const override
HestonSLVFDMModel(Handle< LocalVolTermStructure > localVol, Handle< HestonModel > hestonModel, const Date &endDate, HestonSLVFokkerPlanckFdmParams params, bool logging=false, std::vector< Date > mandatoryDates=std::vector< Date >(), Real mixingFactor=1.0)
const std::list< LogEntry > & logEntries() const
const HestonSLVFokkerPlanckFdmParams params_
ext::shared_ptr< LocalVolTermStructure > leverageFunction_
std::list< LogEntry > logEntries_
const Handle< LocalVolTermStructure > localVol_
const Handle< HestonModel > hestonModel_
const std::vector< Date > mandatoryDates_
ext::shared_ptr< HestonProcess > hestonProcess() const
ext::shared_ptr< LocalVolTermStructure > localVol() const
virtual void calculate() const
Linear-interpolation factory and traits
std::vector< Size > rescaleTimeSteps() const
ext::shared_ptr< Fdm1dMesher > mesher(Time t) const
Real invcdf(Real p, Time t) const override
Matrix used in linear algebra.
template class providing a null value for a given type.
std::pair< iterator, bool > registerWith(const ext::shared_ptr< Observable > &)
Real Time
continuous quantity with 1-year units
Real Volatility
volatility
std::size_t Size
size of a container
Array Pow(const Array &v, Real alpha)
Array Exp(const Array &v)
static FdmSchemeDesc ImplicitEuler()
const Size tMinStepsPerYear
const Real tStepNumberDecay
const Size tMaxStepsPerYear
const Size maxIntegrationIterations
const FdmSchemeDesc schemeDesc
const Real localVolEpsProb
const Size nRannacherTimeSteps
const Size predictionCorretionSteps
const Real leverageFctPropEps
const FdmSquareRootFwdOp::TransformationType trafoType
const FdmHestonGreensFct::Algorithm greensAlgorithm