QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.38
Loading...
Searching...
No Matches
operatorsplittingspreadengine.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) 2024 Klaus Spanderen
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
20#include <ql/math/functional.hpp>
21#include <ql/math/distributions/normaldistribution.hpp>
22#include <ql/pricingengines/basket/operatorsplittingspreadengine.hpp>
23
24#include <utility>
25
26namespace QuantLib {
27
28 OperatorSplittingSpreadEngine::OperatorSplittingSpreadEngine(
29 ext::shared_ptr<GeneralizedBlackScholesProcess> process1,
30 ext::shared_ptr<GeneralizedBlackScholesProcess> process2,
31 Real correlation,
32 Order order)
33 : SpreadBlackScholesVanillaEngine(std::move(process1), std::move(process2), correlation),
34 order_(order) {
35 }
36
37 Real OperatorSplittingSpreadEngine::calculate(
38 Real f1, Real f2, Real k, Option::Type optionType,
39 Real variance1, Real variance2, DiscountFactor df) const {
40
41 const auto callPutParityPrice = [f1, f2, df, k, optionType](Real callPrice) -> Real {
42 if (optionType == Option::Call)
43 return callPrice;
44 else
45 return callPrice - df*(f1-f2-k);
46 };
47
48 const Real vol1 = std::sqrt(variance1);
49 const Real vol2 = std::sqrt(variance2);
50 const Real sig2 = vol2*f2/(f2+k);
51 const Real sig_m = std::sqrt(variance1 +sig2*(sig2 - 2*rho_*vol1));
52
53 const Real d1 = (std::log(f1) - std::log(f2 + k))/sig_m + 0.5*sig_m;
54 const Real d2 = d1 - sig_m;
55
56 const CumulativeNormalDistribution N;
57
58 const Real kirkCallNPV = df*(f1*N(d1) - (f2 + k)*N(d2));
59
60 const Real vs = vol2/(sig_m*sig_m);
61 const Real rs = squared(rho_*vol1 - sig2);
62
63 const Real oPlt = -sig2*sig2*k*df*NormalDistribution()(d2)*vs
64 *( -d2*rs/sig2 - 0.5*vs*sig_m * k / (f2+k)
65 *(rs*d1*d2 + (1-rho_*rho_)*variance1));
66
67 if (order_ == First)
68 return callPutParityPrice(kirkCallNPV + 0.5*oPlt);
69
70 QL_REQUIRE(order_ == Second, "unknown approximation type");
71
72 /*
73 In the original paper, the second-order approximation was computed using numerical differentiation.
74 The following Mathematica scripts calculates the approximation to the n'th order.
75
76 vol2Hat[R2_] := vol2*(R2 - K)/R2
77 volMinusHat[R2_] := Sqrt[vol1^2 - 2*rho*vol1*vol2Hat[R2] + vol2Hat[R2]^2]
78 zeta1[R1_, R2_] := 1/(volMinusHat[R2]*Sqrt[t])*(Log[R1] + volMinusHat[R2]^2*t/2)
79 zeta2[R1_, R2_] := zeta1[R1, R2] - volMinusHat[R2]*Sqrt[t]
80 pLT[R1_, R2_] := Exp[-r*t]*R2*(R1*CDF[NormalDistribution[0, 1], zeta1[R1, R2]]
81 - CDF[NormalDistribution[0, 1], zeta2[R1, R2]])
82 opt[R1_, R2_] := (1/2*vol2Hat[R2]^2*R2^2*D[#, {R2, 2}] + (rho*vol1 - vol2Hat[R2])*vol2Hat[R2]*R1*R2*
83 D[#, R1, R2] - (rho*vol1 - vol2Hat[R2])*vol2Hat[R2]*R1*D[#, R1]) &
84
85 pStrange1[R1_, R2_] := pLT[R1, R2] + (t/2)^1/Factorial[1]*opt[R1, R2][pLT[R1, R2]]
86 pStrange2[R1_, R2_] := pStrange1[R1, R2] + (t/2)^2/Factorial[2]*opt[R1, R2][opt[R1, R2][pLT[R1, R2]]]
87 */
88
89 const Real R2 = f2+k;
90 const Real R1 = f1/R2;
91 const Real vol12 = vol1*vol1;
92 const Real vol22 = vol2*vol2;
93 const Real vol23 = vol22*vol2;
94
95 if (rs < std::pow(QL_EPSILON, 0.625)) {
96 const Real vol24 = vol22*vol22;
97 const Real vol26 = vol22*vol24;
98 const Real k2 = k*k;
99 const Real R22 = R2*R2;
100 const Real R24 = R22*R22;
101 const Real kmR22 = squared(k-R2);
102 const Real lnR1 = std::log(R1);
103
104 const Real ooPlt = -0.0625*(k2*kmR22*vol26*(-8*R22*R24*(7*k2 - 7*k*R2 + R22)*vol12*vol12
105 + kmR22*R24*vol12*(-112*k*R2 + 16*R22 + k2*(124 + 3*vol12))*vol22 - 2*kmR22*kmR22*R22*
106 (-28*k*R2 + 4*R22 + k2*(34 + 3*vol12))*vol24 + 3*k2*kmR22*kmR22*kmR22*vol26 - 4*k*(k -R2)*
107 R24*lnR1*(-4*R22*vol12 + 4*kmR22*vol22 + 3*k*(k - R2)*vol22*lnR1)))
108 /(std::exp(squared(-(R22*vol12) + kmR22*vol22 + 2*R22*lnR1)/(8*R24*vol12 - 8*kmR22*R22*vol22))*
109 M_SQRTPI*M_SQRT2*R22*R24*R2*squared(R22*vol12 - kmR22*vol22)*std::sqrt(vol12 - (kmR22*vol22)/R22));
110
111 return callPutParityPrice(kirkCallNPV + 0.5*oPlt + 0.125*ooPlt);
112 }
113
114 const Real F2 = f2;
115
116 const Real F22 = F2*F2;
117 const Real F23 = F22*F2;
118 const Real F24 = F22*F22;
119
120 const Real iR2 = 1.0/R2;
121 const Real iR22 = iR2*iR2;
122 const Real iR23 = iR22*iR2;
123 const Real iR24 = iR22*iR22;
124 const Real a = vol12 - 2*F2*iR2*rho_*vol1*vol2 + F22*iR22*vol22;
125 const Real a2 = a*a;
126 const Real b = a/2+std::log(R1);
127 const Real b2 = b*b;
128 const Real c = std::sqrt(a);
129 const Real d = b/c;
130 const Real e = rho_*vol1 - F2*iR2*vol2;
131 const Real e2 = e*e;
132 const Real f = d-c;
133 const Real g = -2*iR2*rho_*vol1*vol2 + 2*F2*iR22*rho_*vol1*vol2 + 2*F2*iR22*vol22 - 2*F22*iR23*vol22;
134 const Real h = rho_*rho_;
135 const Real j = 1-h;
136 const Real iat = 1/c;
137 const Real l = b*iat - c;
138 const Real m = f*(1 - (R2*rho_*vol1)/(F2*vol2)) - (e*iR2*k*(d*l + (j*vol12)/(e*e))*vol2)/(2.*c);
139 const Real n = (iat*(1 - (R2*rho_*vol1)/(F2*vol2)))/R1 - (e*iR2*k*((f*iat)/R1 + b/(a*R1))*vol2)/(2.*c);
140 const Real o = df*std::exp(-0.5*f*f);
141 const Real p = d*l + (j*vol12)/(e*e);
142 const Real q = (-2*j*vol12*(-(iR2*vol2) + F2*iR22*vol2))/(e*e*e);
143 const Real s = q - (b2*g)/(2.*a2) - (b*f*g)/(2.*a*c) + (f*g)/(2.*c);
144 const Real u = f*(-((rho_*vol1)/(F2*vol2)) + (R2*rho_*vol1)/(F22*vol2));
145 const Real v = -0.5*(b*g*(1 - (R2*rho_*vol1)/(F2*vol2)))/(a*c);
146 const Real w = (3*g*g)/(4.*a2*c) - (4*iR22*rho_*vol1*vol2 - 4*F2*iR23*rho_*vol1*vol2 +
147 2*iR22*vol22 - 8*F2*iR23*vol22 + 6*F22*iR24*vol22)/(2.*a*c);
148 const Real x = u + v + (e*g*iR2*k*p*vol2)/(4.*a*c) + (e*iR22*k*p*vol2)/(2.*c) -
149 (e*iR2*k*s*vol2)/(2.*c) - (iR2*k*p*vol2*(-(iR2*vol2) + F2*iR22*vol2))/(2.*c);
150 const Real y = (4*iR22 - 4*F2*iR23)*rho_*vol1*vol2 + (2*iR22 - 8*F2*iR23 + 6*F22*iR24)*vol22;
151 const Real z = 4*iR22*rho_*vol1*vol2 - 4*F2*iR23*rho_*vol1*vol2 + 2*iR22*vol22 -
152 8*F2*iR23*vol22 + 6*F22*iR24*vol22;
153
154 const Real ooPlt = (k*o*vol23*(-2*c*b2*e2*e*(-1 + f*f)*F23*F24*g*g*iR22*m*vol23 +
155 2*b2*e2*e2*F23*F24*g*g*iR2*iR22*k*vol22*vol22 + 2*a*b*e2*e*F23*F22*g*iR22*vol2*(-8*
156 e2*F2*iR2*k*vol22 + 7*f*F22*g*m*vol22) - a*c*e2*e*F23*F22*g*iR22*vol2*(4*e*F2*
157 vol2*(-2*b*(-1 + f*f)*m + e*f*iR2*k*vol2) + F22*g*(16*m + e*(2*f + 3*b*iat)*iR2*k*vol2)*vol22) -
158 4*a2*a*c*e2*(e2*F22*vol2*(4*F22*iat*iR22*R2*rho_*vol1 + 8*F23*iR22*n*R1*vol2 -
159 4*F24*3*iR23*n*R1*vol2 - F23*iR22*(4*iat*rho_*vol1 + F22*iR2*k*p*vol23*w))
160 + 4*F23*F22*vol22*vol22*(iR22*(-2*F2*iR2 + 3*F22*iR22)*m + F22*(2*iR2 - 3*F2*iR22)*iR23*m
161 + F22*iR22*(-iR2 + F2*iR22)*x) + 2*e*F22*(2*F24*F2*iR24*n*R1*vol23 +
162 2*f*F2*F22*iR22*rho_*vol1*vol22 - 2*f*F22*iR22*R2*rho_*vol1*vol22 -
163 b*F24*iR22*R2*rho_*vol1*vol22*w - 2*F24*vol2*(iR23*n*R1*vol22 + 4*iR23*m*vol22 - 2*iR22*vol22*x) +
164 F23*(2*iR22*m*vol23 + 6*F22*iR24*m*vol23 + b*F22*iR22*vol23*w - 4*F22*iR23*vol23*x))) +
165 2*a2*c*e2*F23*F22*vol2*(8*F22*g*iR22*(-iR2 + F2*iR22)*m*vol23 +
166 e2*iR22*vol2*(8*F2*g*n*R1 + b*F22*iat*iR2*k*vol22*(y - z)) +
167 4*e*vol22*(4*F2*g*iR22*m + F22*(-4*g*iR23*m + 2*g*iR22*x + iR22*m*z))) +
168 2*a2*a*F22*(-4*e2*e2*e*f*F24*iat*iR24*k*vol23 + 8*e*F2*F24*iR23*(-iR22 + F2*iR23)*j*k*vol12*vol23*vol22 +
169 12*F2*F24*iR23*squared(iR2 - F2*iR22)*j*k*vol12*vol23*vol23 +
170 e2*e2*F2*vol22*(2*F24*iR22*k*vol22*(2*(iR23*p - iR22*s) + b2*iat*iR2*w) + f*(4*F22*iR22*(4*m +
171 F22*iat*iR2*iR22*k*vol22) - 4*F23*(6*iR23*m + iat*iR24*k*vol22 - 2*iR22*x)
172 + F24*iR23*k*vol22*(2*b*w + iat*y))) - 2*e2*e*F22*iR22*(4*f*F22*(iR2 - F2*iR22)*m*vol23 +
173 F22*vol22*(F2*vol2*(2*k*(F2*iR24*p + F2*iR24*p + iR22*s - iR23*(2*p + F2*s))*vol22 + y - z)
174 + R2*rho_*vol1*(-y + z)))) - 2*a2*e2*F23*(2*e2*e*F23*iR22*k*(2*b*iR22 + g*(-1 + f*iat)*iR2)*vol23 +
175 4*b*f*F22*F22*g*iR22*(-iR2 + F2*iR22)*m*vol22*vol22 +
176 2*e2*F22*iR22*vol2*(2*b*F2*iR2*(iR2 - F2*iR22)*k*vol23 + g*(2*R2*rho_*vol1 + 2*F2*(-1 + 3*f*m +
177 b*f*n*R1)*vol2 + F22*k*(-(iR22*p) + iR2*s)*vol23)) + e*vol22*(F2*F22*g*iR22*(g*R2*rho_*vol1 + F2*g*(-1 + f*m)*vol2 +
178 2*F2*iR2*(-iR2 + F2*iR22)*k*p*vol23) + 2*b*(2*F2*F22*g*iR22*rho_*vol1 - 2*F22*g*iR22*R2*rho_*vol1 +
179 4*f*F23*g*iR22*m*vol2 + f*F24*vol2*(-4*g*iR23*m + 2*g*iR22*x +
180 iR22*m*z))))))/(16.*a2*a2*c*e2*F23*M_SQRT2*M_SQRTPI*vol2);
181
182 return callPutParityPrice(kirkCallNPV + 0.5*oPlt + 0.125*ooPlt);
183 }
184}
185
186
QuantLib::CumulativeNormalDistribution
Cumulative normal distribution function.
Definition: normaldistribution.hpp:68
QuantLib::NormalDistribution
Normal distribution function.
Definition: normaldistribution.hpp:44
QuantLib::OperatorSplittingSpreadEngine::OperatorSplittingSpreadEngine
OperatorSplittingSpreadEngine(ext::shared_ptr< GeneralizedBlackScholesProcess > process1, ext::shared_ptr< GeneralizedBlackScholesProcess > process2, Real correlation, Order order=Second)
Definition: operatorsplittingspreadengine.cpp:28
QuantLib::PricingEngine::calculate
virtual void calculate() const =0
QL_REQUIRE
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
Definition: errors.hpp:117
d
Date d
Definition: exchangeratemanager.cpp:32
b
std::function< Real(Real)> b
Definition: extendedornsteinuhlenbeckprocess.cpp:30
QL_EPSILON
#define QL_EPSILON
Definition: qldefines.hpp:178
QuantLib::Real
QL_REAL Real
real number
Definition: types.hpp:50
QuantLib::DiscountFactor
Real DiscountFactor
discount factor between dates
Definition: types.hpp:66
order_
const Size order_
Definition: lsmbasissystem.cpp:51
functional.hpp
functionals and combinators not included in the STL
M_SQRTPI
#define M_SQRTPI
Definition: mathconstants.hpp:74
M_SQRT2
#define M_SQRT2
Definition: mathconstants.hpp:94
QuantLib
Definition: any.hpp:37
QuantLib::squared
T squared(T x)
Definition: functional.hpp:37
std
STL namespace.
normaldistribution.hpp
normal, cumulative and inverse cumulative distributions
operatorsplittingspreadengine.hpp
Analytic operator splitting approximation by Chi-Fai Lo (2015)
q
ext::shared_ptr< YieldTermStructure > q
Definition: perturbativebarrieroptionengine.cpp:1464
v
ext::shared_ptr< BlackVolTermStructure > v
Definition: perturbativebarrieroptionengine.cpp:1487
Generated by Doxygen 1.9.5