QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
zabr.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) 2014 Peter Caspers
5
6 This file is part of QuantLib, a free-software/open-source library
7 for financial quantitative analysts and developers - http://quantlib.org/
8
9 QuantLib is free software: you can redistribute it and/or modify it
10 under the terms of the QuantLib license. You should have received a
11 copy of the license along with this program; if not, please email
12 <quantlib-dev@lists.sf.net>. The license is also available online at
13 <http://quantlib.org/license.shtml>.
14
15 This program is distributed in the hope that it will be useful, but WITHOUT
16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17 FOR A PARTICULAR PURPOSE. See the license for more details.
18*/
19
22#include <ql/errors.hpp>
37
38using std::pow;
39
40namespace QuantLib {
41
42ZabrModel::ZabrModel(const Real expiryTime, const Real forward,
43 const Real alpha, const Real beta, const Real nu,
44 const Real rho, const Real gamma)
45 : expiryTime_(expiryTime), forward_(forward), alpha_(alpha), beta_(beta),
46 nu_(nu * std::pow(alpha_, 1.0 - gamma)), rho_(rho), gamma_(gamma) {
47
49 QL_REQUIRE(gamma >= 0.0 /*&& gamma<=1.0*/,
50 "gamma must be non negative: " << gamma << " not allowed");
51 QL_REQUIRE(forward >= 0.0,
52 "forward must be non negative: " << forward << " not allowed");
53 QL_REQUIRE(expiryTime > 0.0, "expiry time must be positive: "
54 << expiryTime << " not allowed");
55}
56
58 const Real x) const {
59 if (close(strike, forward_))
60 return std::pow(forward_, beta_ - 1.0) * alpha_;
61 else
62 return std::log(forward_ / strike) / x;
63}
64
66 return lognormalVolatility(std::vector<Real>(1, strike))[0];
67}
68
69std::vector<Real> ZabrModel::lognormalVolatility(const std::vector<Real> &strikes) const {
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); });
74 return result;
75}
76
77Real ZabrModel::normalVolatilityHelper(const Real strike, const Real x) const {
78 if (close(strike, forward_))
79 return std::pow(forward_, beta_) * alpha_;
80 else
81 return (forward_ - strike) / x;
82}
83
85 return normalVolatility(std::vector<Real>(1, strike))[0];
86}
87
88std::vector<Real> ZabrModel::normalVolatility(const std::vector<Real> &strikes) const {
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); });
93 return result;
94}
95
97 return alpha_ * std::pow(std::fabs(f), beta_) /
98 F(y(f), std::pow(alpha_, gamma_ - 1.0) *
99 x); // TODO optimize this, y is comoputed together
100 // with x already
101}
102
104 return localVolatility(std::vector<Real>(1, f))[0];
105}
106
107std::vector<Real> ZabrModel::localVolatility(const std::vector<Real> &f) const {
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); });
112 return result;
113}
114
115Real ZabrModel::fdPrice(const Real strike) const {
116 return fdPrice(std::vector<Real>(1, strike))[0];
117}
118
119std::vector<Real> ZabrModel::fdPrice(const std::vector<Real> &strikes) const {
120
121 // TODO check strikes to be increasing
122 // TODO put these parameters somewhere
123 const Real start =
124 std::min(0.00001, strikes.front() * 0.5); // lowest strike for grid
125 const Real end =
126 std::max(0.10, strikes.back() * 1.5); // highest strike for grid
127 const Size size = 500; // grid points
128 const Real density = 0.1; // density for concentrating mesher
129 const Size steps =
130 (Size)std::ceil(expiryTime_ * 24); // number of steps in dimension t
131 const Size dampingSteps = 5; // thereof damping steps
132
133#if defined(__GNUC__) && (__GNUC__ >= 12)
134#pragma GCC diagnostic push
135#pragma GCC diagnostic ignored "-Warray-bounds"
136#endif
137
138 // Layout
139 std::vector<Size> dim(1, size);
140 const ext::shared_ptr<FdmLinearOpLayout> layout(
141 new FdmLinearOpLayout(dim));
142
143#if defined(__GNUC__) && (__GNUC__ >= 12)
144#pragma GCC diagnostic pop
145#endif
146
147 // Mesher
148 const ext::shared_ptr<Fdm1dMesher> m1(new Concentrating1dMesher(
149 start, end, size, std::pair<Real, Real>(forward_, density), true));
150 // const ext::shared_ptr<Fdm1dMesher> m1(new
151 // Uniform1dMesher(start,end,size));
152 // const ext::shared_ptr<Fdm1dMesher> m1a(new
153 // Uniform1dMesher(start,0.03,101));
154 // const ext::shared_ptr<Fdm1dMesher> m1b(new
155 // Uniform1dMesher(0.03,end,100));
156 // const ext::shared_ptr<Fdm1dMesher> m1(new Glued1dMesher(*m1a,*m1b));
157 const std::vector<ext::shared_ptr<Fdm1dMesher> > meshers(1, m1);
158 const ext::shared_ptr<FdmMesher> mesher(
159 new FdmMesherComposite(layout, meshers));
160
161 // Boundary conditions
162 FdmBoundaryConditionSet boundaries;
163
164 // initial values
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);
169 }
170
171 // local vols (TODO how can we avoid these Array / vector copying?)
172 Array k = mesher->locations(0);
173 std::vector<Real> kv(k.size());
174 std::copy(k.begin(), k.end(), kv.begin());
175 std::vector<Real> locVolv = localVolatility(kv);
176 Array locVol(locVolv.size());
177 std::copy(locVolv.begin(), locVolv.end(), locVol.begin());
178
179 // solver
180 ext::shared_ptr<FdmDupire1dOp> map(new FdmDupire1dOp(mesher, locVol));
181 FdmBackwardSolver solver(map, boundaries,
182 ext::shared_ptr<FdmStepConditionComposite>(),
184 solver.rollback(rhs, expiryTime_, 0.0, steps, dampingSteps);
185
186 // interpolate solution
187 ext::shared_ptr<Interpolation> solution(new CubicInterpolation(
188 k.begin(), k.end(), rhs.begin(), CubicInterpolation::Spline, true,
191 // ext::shared_ptr<Interpolation> solution(new
192 // LinearInterpolation(k.begin(),k.end(),rhs.begin()));
193 solution->disableExtrapolation();
194 std::vector<Real> result(strikes.size());
195 std::transform(strikes.begin(), strikes.end(), result.begin(), *solution);
196 return result;
197}
198
199Real ZabrModel::fullFdPrice(const Real strike) const {
200
201 // TODO what are good values here, still experimenting with them
202 Real eps = 0.01;
203 Real scaleFactor = 1.5;
204 Real normInvEps = InverseCumulativeNormal()(1.0 - eps);
205 Real alphaI = alpha_ * std::pow(forward_, beta_ - 1.0);
206 // nu is already standardized within this class ...
207 Real v0 = alpha_ * std::exp(-scaleFactor * normInvEps *
208 std::sqrt(expiryTime_) * nu_);
209 Real v1 = alpha_ *
210 std::exp(scaleFactor * normInvEps * std::sqrt(expiryTime_) * nu_);
211 Real f0 = forward_ * std::exp(-scaleFactor * normInvEps *
212 std::sqrt(expiryTime_) * alphaI);
213 Real f1 = forward_ * std::exp(scaleFactor * normInvEps *
214 std::sqrt(expiryTime_) * alphaI);
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)));
218
219 const Size sizef = 100;
220 const Size sizev = 100;
221 const Size steps = Size(24 * expiryTime_ + 1);
222 const Size dampingSteps = 5;
223 const Real densityf = 0.1;
224 const Real densityv = 0.1;
225
226 QL_REQUIRE(strike >= f0 && strike <= f1,
227 "strike (" << strike << ") must be inside pde grid [" << f0
228 << ";" << f1 << "]");
229
230 // Layout
231 std::vector<Size> dim;
232 dim.push_back(sizef);
233 dim.push_back(sizev);
234 const ext::shared_ptr<FdmLinearOpLayout> layout(
235 new FdmLinearOpLayout(dim));
236
237 // Mesher
238 // two concentrating mesher around f and k to get the mesher for f
239 const Real x0 = std::min(forward_, strike);
240 const Real x1 = std::max(forward_, strike);
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; // common point, so we can spend
244 // one more here
245 const ext::shared_ptr<Fdm1dMesher> mfa(
246 new Concentrating1dMesher(f0, (x0 + x1) / 2.0, sizefa,
247 std::pair<Real, Real>(x0, densityf), true));
248 const ext::shared_ptr<Fdm1dMesher> mfb(
249 new Concentrating1dMesher((x0 + x1) / 2.0, f1, sizefb,
250 std::pair<Real, Real>(x1, densityf), true));
251 const ext::shared_ptr<Fdm1dMesher> mf(new Glued1dMesher(*mfa, *mfb));
252
253 // concentraing mesher around f to get the forward mesher
254 // const ext::shared_ptr<Fdm1dMesher> mf(new Concentrating1dMesher(
255 // f0, f1, sizef, std::pair<Real, Real>(forward_, densityf), true));
256
257 // Volatility mesher
258 const ext::shared_ptr<Fdm1dMesher> mv(new Concentrating1dMesher(
259 v0, v1, sizev, std::pair<Real, Real>(alpha_, densityv), true));
260
261 // uniform meshers
262 // const ext::shared_ptr<Fdm1dMesher> mf(new
263 // Uniform1dMesher(f0,f1,sizef));
264 // const ext::shared_ptr<Fdm1dMesher> mv(new
265 // Uniform1dMesher(v0,v1,sizev));
266
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(
271 new FdmMesherComposite(layout, meshers));
272
273 // initial values
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);
279 // Real v = 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));
285 }
286
287 // Boundary conditions
288 FdmBoundaryConditionSet boundaries;
289
290 ext::shared_ptr<FdmZabrOp> map(
291 new FdmZabrOp(mesher, beta_, nu_, rho_, gamma_));
292 FdmBackwardSolver solver(map, boundaries,
293 ext::shared_ptr<FdmStepConditionComposite>(),
294 FdmSchemeDesc::/*CraigSneyd()*/ Hundsdorfer());
295
296 solver.rollback(rhs, expiryTime_, 0.0, steps, dampingSteps);
297
298 // interpolate solution (this is not necessary when using concentrating
299 // meshers with required point)
300 Matrix result(f_.size(), v_.size());
301 for (Size j = 0; j < v_.size(); ++j)
302 std::copy(rhs.begin() + j * f_.size(),
303 rhs.begin() + (j + 1) * f_.size(), result.row_begin(j));
304 ext::shared_ptr<BicubicSpline> interpolation =
305 ext::make_shared<BicubicSpline>(
306 f_.begin(), f_.end(), v_.begin(), v_.end(), result);
307 interpolation->disableExtrapolation();
308 return (*interpolation)(forward_, alpha_);
309}
310
311Real ZabrModel::x(const Real strike) const {
312 return x(std::vector<Real>(1, strike))[0];
313}
314
315std::vector<Real> ZabrModel::x(const std::vector<Real> &strikes) const {
316
317 QL_REQUIRE(strikes[0] > 0.0 || beta_ < 1.0,
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 << ")");
322
323 AdaptiveRungeKutta<Real> rk(1.0E-8, 1.0E-5,
324 0.0); // TODO move the parameters here as
325 // parameters with default values to
326 // the constructor
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); });
330
331 if (close(gamma_, 1.0)) {
332 for (Size m = 0; m < y.size(); m++) {
333 Real J = std::sqrt(1.0 + nu_ * nu_ * y[m] * y[m] -
334 2.0 * rho_ * nu_ * y[m]);
335 result[y.size() - 1 - m] =
336 std::log((J + nu_ * y[m] - rho_) / (1.0 - rho_)) / nu_;
337 }
338 } else {
339 Size ynz = std::upper_bound(y.begin(), y.end(), 0.0) - y.begin();
340 if (ynz > 0)
341 if (close(y[ynz - 1], 0.0))
342 ynz--;
343 if (ynz == y.size())
344 ynz--;
345
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) {
350 Real u = rk([&](Real _y, Real _u){ return F(_y, _u); },
351 u0, y0, y[m]);
352 result[y.size() - 1 - m] = u * pow(alpha_, 1.0 - gamma_);
353 u0 = u;
354 y0 = y[m];
355 }
356 }
357 }
358
359 return result;
360}
361
362Real ZabrModel::y(const Real strike) const {
363
364 if (close(beta_, 1.0)) {
365 return std::log(forward_ / strike) * std::pow(alpha_, gamma_ - 2.0);
366 } else {
367 return (strike < 0.0
368 ? Real(std::pow(forward_, 1.0 - beta_) +
369 std::pow(-strike, 1.0 - beta_))
370 : Real(std::pow(forward_, 1.0 - beta_) -
371 std::pow(strike, 1.0 - beta_))) *
372 std::pow(alpha_, gamma_ - 2.0) / (1.0 - beta_);
373 }
374}
375
376Real ZabrModel::F(const Real y, const Real u) const {
377 Real A = 1.0 + (gamma_ - 2.0) * (gamma_ - 2.0) * nu_ * nu_ * y * y +
378 2.0 * rho_ * (gamma_ - 2.0) * nu_ * y;
379 Real B = 2.0 * rho_ * (1.0 - gamma_) * nu_ +
380 2.0 * (1.0 - gamma_) * (gamma_ - 2.0) * nu_ * nu_ * y;
381 Real C = (1.0 - gamma_) * (1.0 - gamma_) * nu_ * nu_;
382 return (-B * u + std::sqrt(B * B * u * u - 4.0 * A * (C * u * u - 1.0))) /
383 (2.0 * A);
384}
385}
Runge-Kutta ODE integration.
1-D array used in linear algebra.
Definition: array.hpp:52
const_iterator end() const
Definition: array.hpp:511
Size size() const
dimension of the array
Definition: array.hpp:495
const_iterator begin() const
Definition: array.hpp:503
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.
Definition: matrix.hpp:41
const_row_iterator row_begin(Size i) const
Definition: matrix.hpp:360
Real x(Real strike) const
Definition: zabr.cpp:311
const Real beta_
Definition: zabr.hpp:67
Real normalVolatility(Real strike) const
Definition: zabr.cpp:84
const Real rho_
Definition: zabr.hpp:67
Real lognormalVolatility(Real strike) const
Definition: zabr.cpp:65
const Real alpha_
Definition: zabr.hpp:67
Real localVolatility(Real f) const
Definition: zabr.cpp:103
Real forward() const
Definition: zabr.hpp:57
Real fdPrice(Real strike) const
Definition: zabr.cpp:115
const Real gamma_
Definition: zabr.hpp:68
Real expiryTime() const
Definition: zabr.hpp:58
Real fullFdPrice(Real strike) const
Definition: zabr.cpp:199
Real localVolatilityHelper(Real f, Real x) const
Definition: zabr.cpp:96
Real rho() const
Definition: zabr.hpp:62
Real nu() const
Definition: zabr.hpp:61
Real lognormalVolatilityHelper(Real strike, Real x) const
Definition: zabr.cpp:57
const Real expiryTime_
Definition: zabr.hpp:66
const Real nu_
Definition: zabr.hpp:67
Real gamma() const
Definition: zabr.hpp:63
Real normalVolatilityHelper(Real strike, Real x) const
Definition: zabr.cpp:77
Real beta() const
Definition: zabr.hpp:60
Real alpha() const
Definition: zabr.hpp:59
Real F(Real y, Real u) const
Definition: zabr.cpp:376
Real y(Real strike) const
Definition: zabr.cpp:362
ZabrModel(Real expiryTime, Real forward, Real alpha, Real beta, Real nu, Real rho, Real gamma)
Definition: zabr.cpp:42
const Real forward_
Definition: zabr.hpp:66
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
Definition: errors.hpp:117
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.
QL_REAL Real
real number
Definition: types.hpp:50
std::size_t Size
size of a container
Definition: types.hpp:58
Real v0
Real rho
Definition: any.hpp:35
bool close(const Quantity &m1, const Quantity &m2, Size n)
Definition: quantity.cpp:163
void validateSabrParameters(Real alpha, Real beta, Real nu, Real rho)
Definition: sabr.cpp:145
OperatorTraits< FdmLinearOp >::bc_set FdmBoundaryConditionSet
STL namespace.
normal, cumulative and inverse cumulative distributions
Differential operator traits.
Real beta
Definition: sabr.cpp:200
Real F
Definition: sabr.cpp:200
Real nu
Definition: sabr.cpp:200
Real alpha
Definition: sabr.cpp:200
SABR functions.
static FdmSchemeDesc Douglas()
One-dimensional simple uniform grid mesher.
std::uint64_t x_
ZABR functions Reference: Andreasen, Huge: ZABR - Expansions for the masses, Preliminary Version,...