45 : expiryTime_(expiryTime), forward_(forward), alpha_(
alpha), beta_(
beta),
46 nu_(
nu *
std::pow(alpha_, 1.0 - gamma)), rho_(
rho), gamma_(gamma) {
50 "gamma must be non negative: " <<
gamma <<
" not allowed");
52 "forward must be non negative: " <<
forward <<
" not allowed");
70 std::vector<Real>
x_ =
x(strikes);
71 std::vector<Real> result(strikes.size());
72 std::transform(strikes.begin(), strikes.end(),
x_.begin(), result.begin(),
73 [&](
Real _k,
Real _x) { return lognormalVolatilityHelper(_k, _x); });
89 std::vector<Real>
x_ =
x(strikes);
90 std::vector<Real> result(strikes.size());
91 std::transform(strikes.begin(), strikes.end(),
x_.begin(), result.begin(),
92 [&](
Real _k,
Real _x) { return normalVolatilityHelper(_k, _x); });
108 std::vector<Real>
x_ =
x(
f);
109 std::vector<Real> result(
f.size());
110 std::transform(
f.begin(),
f.end(),
x_.begin(), result.begin(),
111 [&](
Real _f,
Real _x) { return localVolatilityHelper(_f, _x); });
116 return fdPrice(std::vector<Real>(1, strike))[0];
124 std::min(0.00001, strikes.front() * 0.5);
126 std::max(0.10, strikes.back() * 1.5);
127 const Size size = 500;
128 const Real density = 0.1;
131 const Size dampingSteps = 5;
133#if defined(__GNUC__) && (__GNUC__ >= 12)
134#pragma GCC diagnostic push
135#pragma GCC diagnostic ignored "-Warray-bounds"
139 std::vector<Size> dim(1, size);
140 const ext::shared_ptr<FdmLinearOpLayout> layout(
143#if defined(__GNUC__) && (__GNUC__ >= 12)
144#pragma GCC diagnostic pop
149 start, end, size, std::pair<Real, Real>(
forward_, density),
true));
157 const std::vector<ext::shared_ptr<Fdm1dMesher> > meshers(1, m1);
158 const ext::shared_ptr<FdmMesher> mesher(
165 Array rhs(mesher->layout()->size());
166 for (
const auto& iter : *layout) {
167 Real k = mesher->location(iter, 0);
168 rhs[iter.index()] = std::max(
forward_ - k, 0.0);
172 Array k = mesher->locations(0);
173 std::vector<Real> kv(k.
size());
174 std::copy(k.
begin(), k.
end(), kv.begin());
176 Array locVol(locVolv.size());
177 std::copy(locVolv.begin(), locVolv.end(), locVol.
begin());
180 ext::shared_ptr<FdmDupire1dOp> map(
new FdmDupire1dOp(mesher, locVol));
182 ext::shared_ptr<FdmStepConditionComposite>(),
193 solution->disableExtrapolation();
194 std::vector<Real> result(strikes.size());
195 std::transform(strikes.begin(), strikes.end(), result.begin(), *solution);
203 Real scaleFactor = 1.5;
210 std::exp(scaleFactor * normInvEps * std::sqrt(
expiryTime_) *
nu_);
211 Real f0 =
forward_ * std::exp(-scaleFactor * normInvEps *
213 Real f1 =
forward_ * std::exp(scaleFactor * normInvEps *
215 v1 = std::min(v1, 2.0);
216 f0 = std::min(strike / 2.0, f0);
217 f1 = std::max(strike * 1.5, std::min(f1, std::max(2.0, strike * 1.5)));
219 const Size sizef = 100;
220 const Size sizev = 100;
222 const Size dampingSteps = 5;
223 const Real densityf = 0.1;
224 const Real densityv = 0.1;
227 "strike (" << strike <<
") must be inside pde grid [" << f0
228 <<
";" << f1 <<
"]");
231 std::vector<Size> dim;
232 dim.push_back(sizef);
233 dim.push_back(sizev);
234 const ext::shared_ptr<FdmLinearOpLayout> layout(
241 const Size sizefa = std::max<Size>(
242 4, (
Size)std::ceil(((x0 + x1) / 2.0 - f0) / (f1 - f0) * (
Real)sizef));
243 const Size sizefb = sizef - sizefa + 1;
245 const ext::shared_ptr<Fdm1dMesher> mfa(
247 std::pair<Real, Real>(x0, densityf),
true));
248 const ext::shared_ptr<Fdm1dMesher> mfb(
250 std::pair<Real, Real>(x1, densityf),
true));
251 const ext::shared_ptr<Fdm1dMesher> mf(
new Glued1dMesher(*mfa, *mfb));
259 v0, v1, sizev, std::pair<Real, Real>(
alpha_, densityv),
true));
267 std::vector<ext::shared_ptr<Fdm1dMesher> > meshers;
268 meshers.push_back(mf);
269 meshers.push_back(mv);
270 const ext::shared_ptr<FdmMesher> mesher(
274 Array rhs(mesher->layout()->size());
275 std::vector<Real>
f_;
276 std::vector<Real> v_;
277 for (
const auto& iter : *layout) {
278 Real f = mesher->location(iter, 0);
280 rhs[iter.index()] = std::max(
f - strike, 0.0);
281 if (iter.coordinates()[1] == 0U)
282 f_.push_back(mesher->location(iter, 0));
283 if (iter.coordinates()[0] == 0U)
284 v_.push_back(mesher->location(iter, 1));
290 ext::shared_ptr<FdmZabrOp> map(
293 ext::shared_ptr<FdmStepConditionComposite>(),
300 Matrix result(
f_.size(), v_.size());
301 for (
Size j = 0; j < v_.size(); ++j)
302 std::copy(rhs.
begin() + j *
f_.size(),
304 ext::shared_ptr<BicubicSpline> interpolation =
305 ext::make_shared<BicubicSpline>(
306 f_.begin(),
f_.end(), v_.begin(), v_.end(), result);
307 interpolation->disableExtrapolation();
312 return x(std::vector<Real>(1, strike))[0];
315std::vector<Real>
ZabrModel::x(
const std::vector<Real> &strikes)
const {
318 "strikes must be positive (" << strikes[0] <<
") if beta = 1");
319 for (
auto i = strikes.begin() + 1; i != strikes.end(); ++i)
320 QL_REQUIRE(*i > *(i - 1),
"strikes must be strictly ascending ("
321 << *(i - 1) <<
"," << *i <<
")");
327 std::vector<Real>
y(strikes.size()), result(strikes.size());
328 std::transform(strikes.rbegin(), strikes.rend(),
y.begin(),
329 [&](
Real _k) { return this->y(_k); });
332 for (
Size m = 0; m <
y.size(); m++) {
335 result[
y.size() - 1 - m] =
339 Size ynz = std::upper_bound(
y.begin(),
y.end(), 0.0) -
y.begin();
341 if (
close(
y[ynz - 1], 0.0))
346 for (
int dir = 1; dir >= -1; dir -= 2) {
347 Real y0 = 0.0, u0 = 0.0;
348 for (
int m = ynz + (dir == -1 ? -1 : 0);
349 dir == -1 ? m >= 0 : m < (int)
y.size(); m += dir) {
369 std::pow(-strike, 1.0 -
beta_))
371 std::pow(strike, 1.0 -
beta_))) *
382 return (-B * u + std::sqrt(B * B * u * u - 4.0 * A * (C * u * u - 1.0))) /
Runge-Kutta ODE integration.
1-D array used in linear algebra.
const_iterator end() const
Size size() const
dimension of the array
const_iterator begin() const
Cubic interpolation between discrete points.
@ SecondDerivative
Match value of second derivative at end.
void rollback(array_type &a, Time from, Time to, Size steps, Size dampingSteps)
Inverse cumulative normal distribution function.
Matrix used in linear algebra.
const_row_iterator row_begin(Size i) const
Real x(Real strike) const
Real normalVolatility(Real strike) const
Real lognormalVolatility(Real strike) const
Real localVolatility(Real f) const
Real fdPrice(Real strike) const
Real fullFdPrice(Real strike) const
Real localVolatilityHelper(Real f, Real x) const
Real lognormalVolatilityHelper(Real strike, Real x) const
Real normalVolatilityHelper(Real strike, Real x) const
Real F(Real y, Real u) const
Real y(Real strike) const
ZabrModel(Real expiryTime, Real forward, Real alpha, Real beta, Real nu, Real rho, Real gamma)
floating-point comparisons
One-dimensional grid mesher concentrating around critical points.
ext::function< Real(Real)> f_
Classes and functions for error handling.
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
One-dimensional simple FDM mesher object working on an index.
Dirichlet boundary conditions for differential operators.
Dupire local volatility pricing operator Note that time is reversed in order to make backward solvers...
memory layout of a fdm linear operator
FdmMesher which is a composite of Fdm1dMesher.
Zabr linear pricing operator.
One-dimensional grid mesher combining two existing ones.
std::size_t Size
size of a container
bool close(const Quantity &m1, const Quantity &m2, Size n)
void validateSabrParameters(Real alpha, Real beta, Real nu, Real rho)
OperatorTraits< FdmLinearOp >::bc_set FdmBoundaryConditionSet
normal, cumulative and inverse cumulative distributions
Differential operator traits.
static FdmSchemeDesc Douglas()
ZABR functions Reference: Andreasen, Huge: ZABR - Expansions for the masses, Preliminary Version,...