24#include <boost/math/special_functions/gamma.hpp>
25#include <boost/functional/hash.hpp>
29class NoArbSabrModel::integrand {
30 const NoArbSabrModel* model;
33 integrand(
const NoArbSabrModel* model,
Real strike)
34 : model(model), strike(strike) {}
36 return std::max(
f - strike, 0.0) * model->p(
f);
40class NoArbSabrModel::p_integrand {
41 const NoArbSabrModel* model;
43 explicit p_integrand(
const NoArbSabrModel* model)
53 : expiryTime_(expiryTime), externalForward_(forward), alpha_(
alpha),
55 numericalForward_(forward) {
58 "expiryTime (" <<
expiryTime <<
") out of bounds");
61 "beta (" <<
beta <<
") out of bounds");
65 "sigmaI = alpha*forward^(beta-1.0) ("
66 << sigmaI <<
") out of bounds, alpha=" <<
alpha
69 "nu (" <<
nu <<
") out of bounds");
71 "rho (" <<
rho <<
") out of bounds");
93 ext::make_shared<GaussLobattoIntegral>(
120 ((*integrator_)(integrand(
this, strike),
121 strike, std::max(
fmax_, 2.0 * strike)) /
131 * ((*integrator_)(p_integrand(
this),
132 strike, std::max(
fmax_, 2.0 * strike)) /
166 Real gamma = 1.0 / (2.0 * (1.0 -
beta_));
169 (
nu_ * zf * std::log(zf * Jz / zF) +
170 (1 +
rho_ *
nu_ * zf) / sqrtOmR *
171 (std::atan((
nu_ * z -
rho_) / sqrtOmR) +
172 std::atan(
rho_ / sqrtOmR)));
176 std::pow(zf, 1.0 - gamma) * std::pow(zF, gamma) *
189 : forward_(forward), expiryTime_(expiryTime), alpha_(
alpha), beta_(
beta),
190 nu_(
nu), rho_(
rho), gamma_(1.0 / (2.0 * (1.0 - beta_))) {
195 0.25, 0.5, 0.75, 1.0, 1.25, 1.5, 1.75, 2.0, 2.25, 2.5, 2.75, 3.0,
196 3.25, 3.5, 3.75, 4.0, 4.25, 4.5, 4.75, 5.0, 5.25, 5.5, 5.75, 6.0, 6.25,
197 6.5, 6.75, 7.0, 7.25, 7.5, 7.75, 8.0, 8.25, 8.5, 8.75, 9.0, 9.25, 9.5,
198 9.75, 10.0, 10.25, 10.5, 10.75, 11.0, 11.25, 11.5, 11.75, 12.0, 12.25,
199 12.5, 12.75, 13.0, 13.25, 13.5, 13.75, 14.0, 14.25, 14.5, 14.75, 15.0,
200 15.25, 15.5, 15.75, 16.0, 16.25, 16.5, 16.75, 17.0, 17.25, 17.5, 17.75,
201 18.0, 18.25, 18.5, 18.75, 19.0, 19.25, 19.5, 19.75, 20.0, 20.25, 20.5,
202 20.75, 21.0, 21.25, 21.5, 21.75, 22.0, 22.25, 22.5, 22.75, 23.0, 23.25,
203 23.5, 23.75, 24.0, 24.25, 24.5, 24.75, 25.0, 25.25, 25.5, 25.75, 26.0,
204 26.25, 26.5, 26.75, 27.0, 27.25, 27.5, 27.75, 28.0, 28.25, 28.5, 28.75,
205 29.0, 29.25, 29.5, 29.75, 30.0
209 1.0, 0.8, 0.7, 0.6, 0.5, 0.45, 0.4, 0.35, 0.3, 0.27, 0.24, 0.21,
210 0.18, 0.15, 0.125, 0.1, 0.075, 0.05
213 rhoG_ = { 0.75, 0.50, 0.25, 0.00, -0.25, -0.50, -0.75 };
215 nuG_ = { 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8 };
217 betaG_ = { 0.01, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9 };
227 if (tauInd ==
tauG_.size())
232 expiryTimeTmp =
tauG_.front();
234 Real tauL = (expiryTimeTmp -
tauG_[tauInd - 1]) /
252 if (rhoInd ==
rhoG_.size()) {
260 if (nuInd ==
nuG_.size())
262 Real tmpNuG = nuInd > 0 ?
nuG_[nuInd - 1] : 0.0;
263 Real nuL = (
nu_ - tmpNuG) / (
nuG_[nuInd] - tmpNuG);
269 if (betaInd ==
betaG_.size())
272 tmpBetaG =
betaG_[betaInd];
277 for (
int iTau = -1; iTau <= 0; ++iTau) {
278 for (
int iSigma = -1; iSigma <= 0; ++iSigma) {
279 for (
int iRho = -1; iRho <= 0; ++iRho) {
280 for (
int iNu = -1; iNu <= 0; ++iNu) {
281 for (
int iBeta = -1; iBeta <= 0; ++iBeta) {
283 if (iNu == -1 && nuInd == 0) {
289 if (iBeta == 0 && betaInd ==
betaG_.size()) {
293 int ind = (tauInd + iTau +
294 (sigmaIInd + iSigma +
296 (nuInd + iNu + ((betaInd + iBeta) *
302 "absorption matrix index ("
303 << ind <<
") invalid");
308 phiRes += phiTmp * (iTau == -1 ? (1.0 - tauL) : tauL) *
309 (iSigma == -1 ? (1.0 - sigmaIL) : sigmaIL) *
310 (iRho == -1 ? (1.0 - rhoL) : rhoL) *
311 (iNu == -1 ? (1.0 - nuL) : nuL) *
312 (iBeta == -1 ? (1.0 - betaL) : betaL);
NoArbSabrModel(Real expiryTime, Real forward, Real alpha, Real beta, Real nu, Real rho)
Real digitalOptionPrice(Real strike) const
Real numericalIntegralOverP_
Real optionPrice(Real strike) const
const Real externalForward_
ext::shared_ptr< GaussLobattoIntegral > integrator_
Real forwardError(Real forward) const
std::vector< Real > sigmaIG_
D0Interpolator(Real forward, Real expiryTime, Real alpha, Real beta, Real nu, Real rho)
std::vector< Real > betaG_
std::vector< Real > tauG_
std::vector< Real > rhoG_
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
ext::function< Real(Real)> b
#define QL_MIN_POSITIVE_REAL
std::size_t Size
size of a container
modified Bessel functions of first and second kind
const Real forward_search_step
const Real density_lower_bound
const Real phiByTau_cutoff
const Real expiryTime_max
const Real forward_accuracy
const Size i_max_iterations
const Real density_threshold
const unsigned long sabrabsprob[1209600]
Real modifiedBesselFunction_i_exponentiallyWeighted(Real nu, Real x)