29#include <boost/accumulators/accumulators.hpp>
30#include <boost/accumulators/statistics/mean.hpp>
31#include <boost/accumulators/statistics/stats.hpp>
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()) {}
51 const ext::shared_ptr<HestonProcess> & process,
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);
60 std::multiset<std::pair<Real, Real> > grid;
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();
67 *(1-std::exp(-process->kappa()*
t))/(4*process->kappa());
69 const Real qMin = 0.0;
70 const Real qMax = std::max(process->v0(),
72 df, ncp, 100, 1e-8)(1-epsilon));
74 const Real minVStep=(qMax-qMin)/(50*
size);
78 grid.insert(std::pair<Real, Real>(qMin, epsilon));
81 ps = (1 - epsilon - p)/(
size-i);
84 df, ncp, 100, 1e-8)(p);
86 const Real vx = std::max(vTmp+minVStep, tmp);
89 grid.insert(std::pair<Real, Real>(vx, p));
93 "something wrong with the grid size");
95 const std::vector<std::pair<Real, Real> > tp(grid.begin(), grid.end());
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);
106 catch (
const Error&) {
108 const Real vol = mixedSigma*
109 std::sqrt(process->theta()/(2*process->kappa()));
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));
117 pGrid[i] = i/(
size-1.0);
118 vGrid[i] = lowerBound + i*(upperBound-lowerBound)/(
size-1.0);
122 Real skewHint = ((process->kappa() != 0.0)
123 ?
Real(std::max(1.0, mixedSigma/process->kappa())) : 1.0);
125 std::sort(pGrid.begin(), pGrid.end());
127 interpolated_volatility(pGrid, vGrid),
128 pGrid.front(), pGrid.back())*std::pow(skewHint, 1.5);
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))
140 std::copy(vGrid.begin(), vGrid.end(),
locations_.begin());
151 const ext::shared_ptr<HestonProcess>& process,
152 const ext::shared_ptr<LocalVolTermStructure>& leverageFct,
158 size, process, maturity, tAvgSteps, epsilon, mixingFactor);
168 if (leverageFct !=
nullptr) {
169 typedef boost::accumulators::accumulator_set<
170 Real, boost::accumulators::stats<
171 boost::accumulators::tag::mean> >
176 const Real s0 = process->s0()->value();
178 acc(leverageFct->localVol(0, s0,
true));
183 for (
Size l=1; l <= tAvgSteps; ++l) {
184 const Real t = (maturity*l)/tAvgSteps;
187 const Real fwd = s0*qTS->discount(
t)/rTS->discount(
t);
189 const Size sAvgSteps = 50;
191 std::vector<Real> u(sAvgSteps), sig(sAvgSteps);
193 for (
Size i=0; i < sAvgSteps; ++i) {
194 u[i] = epsilon + ((1-2*epsilon)/(sAvgSteps-1))*i;
197 const Real gf = x*vol*std::sqrt(
t);
198 const Real f = fwd*std::exp(gf);
200 sig[i] =
squared(leverageFct->localVol(
t,
f,
true));
203 const Real leverageAvg =
205 interpolated_volatility(u, sig), u.front(), u.back())
Chi-square (central and non-central) distributions.
std::vector< Real > locations_
Real dminus(Size index) const
Real location(Size index) const
std::vector< Real > dplus_
std::vector< Real > dminus_
Real dplus(Size index) const
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)
Real volaEstimate() const
Integral of a one-dimensional function.
Shared handle to an observable.
Inverse cumulative normal distribution function.
template class providing a null value for a given type.
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
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
std::size_t Size
size of a container
linear interpolation between discrete points
Local volatility term structure base class.
functionals and combinators not included in the STL
normal, cumulative and inverse cumulative distributions