25#include <ql/pricingengines/vanilla/analyticgjrgarchengine.hpp>
26#include <ql/math/distributions/normaldistribution.hpp>
27#include <ql/instruments/payoffs.hpp>
37 const ext::shared_ptr<GJRGARCHModel>& model)
45 "not an European option");
48 ext::shared_ptr<StrikedTypePayoff> payoff =
49 ext::dynamic_pointer_cast<StrikedTypePayoff>(
arguments_.payoff);
50 QL_REQUIRE(payoff,
"non-striked payoff given");
52 const ext::shared_ptr<GJRGARCHProcess>& process =
model_->process();
54 const Rate riskFreeDiscount = process->riskFreeRate()->discount(
56 const Rate dividendDiscount = process->dividendYield()->discount(
58 const Real spotPrice = process->s0()->value();
59 QL_REQUIRE(spotPrice > 0.0,
"negative or null underlying given");
60 const Real strikePrice = payoff->strike();
62 Size T =
Size(std::lround(process->daysPerYear()*term));
63 Real r = -std::log(riskFreeDiscount/dividendDiscount)/(process->daysPerYear()*term);
64 Real h1 = process->v0();
65 Real b0 = process->omega();
66 Real b2 = process->alpha();
67 Real b1 = process->beta();
68 Real b3 = process->gamma();
69 Real la = process->lambda();
71 Real n = std::exp(-la*la/2)/(M_SQRTPI*M_SQRT2);
72 const Real s = spotPrice;
73 const Real x = strikePrice;
74 Real m1, m2, m3, v1, v2, v3, z1, z2, x1;
75 Real ex, ex2, ex3, ex4;
76 Real sEh = 0.0, sEh2 = 0.0, sEhh = 0.0, sEh1_2eh = 0.0;
77 Real sEhhh = 0.0, sEh2h = 0.0, sEhh2 = 0.0, sEh3 = 0.0;
78 Real sEh1_2eh2 = 0.0, sEh3_2eh = 0.0, sEh1_2ehh = 0.0, sEhh1_2eh = 0.0;
79 Real sEhe2h = 0.0, sEh1_2eh1_2eh = 0.0;
82 Real ST1, ST2, ST3, ST4;
85 Real stdev, sigma, k3, k4;
86 Real d, del, d_, C, A3, A4, Capp;
87 bool constants_match =
false;
91 m1 = b1 + (b2+b3*N)*(1+la*la) + b3*la*n;
92 m2 = b1*b1 + b2*b2*(pow(la,4)+6*la*la+3)
93 + (b3*b3+2*b2*b3)*( pow(la,4)*N
94 +pow(la,3)*n+6*la*la*N+5*la*n+3*N)
95 + 2*b1*b2*(1+la*la) + 2*b3*b1*(la*la*N+la*n+N);
97 + (3*b3*b3*b1+6*b1*b2*b3)*(pow(la,3)*n+5*la*n+3*N
98 +pow(la,4)*N+6*la*la*N)
99 + pow(b2,3)*(15+pow(la,6)+15*pow(la,4)+45*la*la)
100 + (pow(b3,3)+3*b2*b2*b3+3*b3*b3*b2)
101 *(pow(la,5)*n+14*pow(la,3)*n+33*la*n+15*N
102 +15*pow(la,4)*N+45*la*la*N+pow(la,6)*N)
103 + 3*b1*b1*b2*(1+la*la) + 3*b1*b1*b3*(la*n+N+la*la*N)
104 + 3*b1*b2*b2*(3+pow(la,4)+6*la*la);
105 v1 = -2*b2*la - 2*b3*(n+la*N);
106 v2 = -4*b2*b2*(3*la+pow(la,3))
107 - (4*b3*b3+8*b2*b3)*(la*la*n+2*n+pow(la,3)*N+3*la*N)
108 - 4*b1*b2*la - 4*b3*b1*(n+la*N);
109 v3 = -12*b3*b1*(b3+2*b2)*(la*la*n+2*n+pow(la,3)*N+3*la*N)
110 - 6*pow(b2,3)*la*(15+pow(la,4)+10*la*la)
111 - 6*b3*(b3*b3+3*b2*b2+3*b3*b2)
112 *(9*la*la*n+8*n+15*la*N+pow(la,4)*n+pow(la,5)*N
114 - 6*b1*b1*b2*la - 6*b3*b1*b1*(n+la*N)
115 - 12*b2*b2*b1*(3*la+std::pow(la,3));
116 z1 = b1 + b2*(3+la*la) + b3*(la*n+3*N+la*la*N);
117 z2 = b1*b1 + b2*b2*(15+pow(la,4)+18*la*la)
118 + (b3*b3+2*b2*b3)*(pow(la,3)*n+17*la*n+15*N
119 +pow(la,4)*N+18*la*la*N)
120 + 2*b1*b2*(3+la*la) + 2*b3*b1*(la*n+3*N+la*la*N);
121 x1 = -6*b2*la - 2*b3*(4*n+3*la*N);
130 constants_match =
true;
134 if (!
init_ || !constants_match || b0 !=
b0_ || h1 !=
h1_ || T !=
T_) {
140 std::unique_ptr<Real[]> m1ai(
new Real[T]);
141 std::unique_ptr<Real[]> m2ai(
new Real[T]);
142 std::unique_ptr<Real[]> m3ai(
new Real[T]);
143 m1ai[0] = m2ai[0] = m3ai[0] = 1.0;
144 for (i=1; i < T; ++i) {
145 m1ai[i] = m1ai[i-1]*m1;
146 m2ai[i] = m2ai[i-1]*m2;
147 m3ai[i] = m3ai[i-1]*m3;
150 for (i = 0; i < T; ++i) {
155 Real m1im2i = m1i-m2i, m1im3i = m1i-m3i, m2im3i = m2i-m3i;
156 Real Eh = b0*(1-m1i)/(1-m1) + m1i*h1;
157 Real Eh2 = b0*b0*((1+m1)*(1-m2i)/(1-m2)
158 - 2*m1*m1im2i/(m1-m2))/(1-m1)
159 + 2*b0*m1*m1im2i*h1/(m1-m2)
161 Real Eh3 = pow(b0,3)*(
163 + 3*m2*((1-m3i)/(1-m3)-m2im3i/(m2-m3))/(1-m2)
164 + 3*m1*((1-m3i)/(1-m3)-m1im3i/(m1-m3))/(1-m1)
166 ((1-m3i)/(1-m3)-m2im3i/(m2-m3))/(1-m2)
167 + (m2im3i/(m2-m3)-m1im3i/(m1-m3))/(m1-m2)
169 + 3*b0*b0*m1*h1*(m1im3i/(m1-m3)
170 +2*m2*(m1im3i/(m1-m3)-m2im3i/(m2-m3))/(m1-m2))
171 + 3*b0*m2*h1*h1*m2im3i/(m2-m3)
173 Real Eh3_2 = .375*std::pow(Eh,-0.5)*Eh2+.625*std::pow(Eh,1.5);
174 Real Eh5_2 = 1.875*std::pow(Eh,0.5)*Eh2-.875*std::pow(Eh,2.5);
178 for (j = 0; j < T-i-1; ++j) {
179 Real Ehh = b0*Eh*(1-m1ai[j+1])/(1-m1)+ Eh2*m1ai[j+1];
180 Real Ehh2 = b0*b0*Eh*((1+m1)*(1-m2ai[j+1])/(1-m2)
182 -m2ai[j+1])/(m1-m2))/(1-m1)
183 + 2*b0*m1*Eh2*(m1ai[j+1]-m2ai[j+1])/(m1-m2)
185 Real Eh2h = b0*Eh2*(1-m1ai[j+1])/(1-m1)
187 Real Eh1_2eh = v1*m1ai[j]*Eh3_2;
188 Real Eh1_2eh2 = 2*b0*v1*(m1ai[j+1]
189 -m2ai[j+1])*Eh3_2/(m1-m2)
191 Real Ehij = b0*(1-m1ai[i+j+1])/(1-m1)
193 Real Ehh3_2 = 0.375*Ehh2/std::sqrt(Ehij)
194 + 0.75*std::sqrt(Ehij)*Ehh
195 - 0.125*std::pow(Ehij,1.5)*Eh;
196 Real Eh3_2eh = v1*m1ai[j]*Eh5_2;
197 Real Eh3_2e3h = x1*m1ai[j]*Eh5_2;
198 Real Eh1_2eh3_2 = 0.375*Eh1_2eh2/std::sqrt(Ehij)
199 + 0.75*std::sqrt(Ehij)*Eh1_2eh;
204 sEh1_2eh2 += Eh1_2eh2;
206 sEhe2h += b0*Eh*(1-m1ai[j+1])/(1-m1)
208 sEh3_2e3h += Eh3_2e3h;
209 for (k = 0; k < T-i-j-2; ++k) {
210 Real Ehhh = b0*Ehh*(1-m1ai[k+1])/(1-m1)
212 Real Eh1_2ehh = b0*Eh1_2eh*(1-m1ai[k+1])/(1-m1)
213 + m1ai[k+1]*Eh1_2eh2;
215 sEh1_2ehh += Eh1_2ehh;
216 sEhh1_2eh += v1*m1ai[k]*Ehh3_2;
217 sEh1_2eh1_2eh += v1*m1ai[k]*Eh1_2eh3_2;
226 ex2 = T*T*r*r - T*r*sEh + 0.25*SD1 + SD2 - SD3;
227 ST1 = 6*sEhhh + (3*sEhh2 + (3*sEh2h + sEh3));
229 ST3 = 2*sEhh1_2eh + (2*sEh1_2ehh + (2*sEh3_2eh + sEh1_2eh2));
230 ST4 = sEhe2h + (sEhh + (sEh2 + 2*sEh1_2eh1_2eh));
231 ex3 = pow(T*r,3) - 1.5*T*T*r*r*sEh
232 + 3*T*r*(SD1/4+SD2-SD3) + (ST2-ST1/8+3*ST3/4-3*ST4/2);
233 SQ2 = 6*sEhe2h + (12*sEh1_2eh1_2eh + 3*sEh2);
234 SQ4 = 2*sEhhh + 2*sEhh2;
235 SQ5 = 3*sEhh1_2eh + 3*sEh1_2ehh + 3*sEh3_2eh
236 + 3*sEh1_2eh2 + sEh3_2e3h;
237 ex4 = pow(T*r,4) - 2*pow(T*r,3)*sEh
238 + 6*T*T*r*r*(SD1/4+SD2-SD3) + T*r*(4*ST2-ST1/2+3*ST3-6*ST4)
239 + (SQ2+3*SQ4/2-2*SQ5);
244 k3 = ex3 - 3*sigma*ex - ex*ex*ex;
246 k4 = ex4 + 6*ex*ex*ex2 - 3*ex*ex*ex*ex - 4*ex*ex3;
247 k3 /= std::pow(sigma,1.5);
257 stdev = std::sqrt(sigma);
258 del = (ex - r*T + sigma/2)/stdev;
259 d = (std::log(s/x) + (r*T+sigma/2))/stdev;
263 A3 = s*std::exp(del*stdev)*stdev*((2*stdev-d_)
264 *std::exp(-d_*d_/2)/std::sqrt(2*M_PI)
266 A4 = s*std::exp(del*stdev)*stdev*(
267 (d_*d_-1-3*stdev*(d_-stdev))*exp(-d_*d_/2)/std::sqrt(2*M_PI)
269 Capp = C + k3*A3 + (k4-3)*A4;
272 switch (payoff->optionType()) {
277 results_.value = Capp+strikePrice*riskFreeDiscount/dividendDiscount
281 QL_FAIL(
"unknown option type");
void calculate() const override
AnalyticGJRGARCHEngine(const ext::shared_ptr< GJRGARCHModel > &model)
Cumulative normal distribution function.
GJR-GARCH model for the stochastic volatility of an asset.
Base class for some pricing engine on a particular model.
Handle< GJRGARCHModel > model_
Vanilla option (no discrete dividends, no barriers) on a single asset.
std::size_t Size
size of a container