1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2008 Yee Man Chan
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/*! \file gjrgarchmodel.hpp
21 \brief analytical approximation pricing engine for a GJR-GARCH option
22 based on Edgeworth expansion
23*/
24
25#include <ql/pricingengines/vanilla/analyticgjrgarchengine.hpp>
26#include <ql/math/distributions/normaldistribution.hpp>
27#include <ql/instruments/payoffs.hpp>
28#include <cmath>
29
30using std::exp;
31using std::pow;
32
33namespace QuantLib {
34
35
36 AnalyticGJRGARCHEngine::AnalyticGJRGARCHEngine(
37 const ext::shared_ptr<GJRGARCHModel>& model)
38 : GenericModelEngine<GJRGARCHModel,
39 VanillaOption::arguments,
40 VanillaOption::results>(model) {init_ = false;}
41
42 void AnalyticGJRGARCHEngine::calculate() const {
43 // this is a european option pricer
44 QL_REQUIRE(arguments_.exercise->type() == Exercise::European,
45 "not an European option");
46
47 // plain vanilla
48 ext::shared_ptr<StrikedTypePayoff> payoff =
49 ext::dynamic_pointer_cast<StrikedTypePayoff>(r: arguments_.payoff);
50 QL_REQUIRE(payoff, "non-striked payoff given");
51
52 const ext::shared_ptr<GJRGARCHProcess>& process = model_->process();
53
54 const Rate riskFreeDiscount = process->riskFreeRate()->discount(
55 d: arguments_.exercise->lastDate());
56 const Rate dividendDiscount = process->dividendYield()->discount(
57 d: arguments_.exercise->lastDate());
58 const Real spotPrice = process->s0()->value();
59 QL_REQUIRE(spotPrice > 0.0, "negative or null underlying given");
60 const Real strikePrice = payoff->strike();
61 const Real term = process->time(arguments_.exercise->lastDate());
62 Size T = Size(std::lround(x: process->daysPerYear()*term));
63 Real r = -std::log(x: 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();
70 Real N = CumulativeNormalDistribution()(la);
71 Real n = std::exp(x: -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;
80 Real sEh3_2e3h = 0.0;
81 Real SD1, SD2, SD3;
82 Real ST1, ST2, ST3, ST4;
83 Real SQ2, SQ4, SQ5;
84 Size i, j, k;
85 Real stdev, sigma, k3, k4;
86 Real d, del, d_, C, A3, A4, Capp;
87 bool constants_match = false;
88
89 if (!init_ || b1 != b1_ || b2 != b2_ || b3 != b3_ || la != la_) {
90 // compute the useful coefficients
91 m1 = b1 + (b2+b3*N)*(1+la*la) + b3*la*n; // ok
92 m2 = b1*b1 + b2*b2*(pow(x: la,y: 4)+6*la*la+3)
93 + (b3*b3+2*b2*b3)*( pow(x: la,y: 4)*N
94 +pow(x: la,y: 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); // ok
96 m3 = pow(x: b1,y: 3)
97 + (3*b3*b3*b1+6*b1*b2*b3)*(pow(x: la,y: 3)*n+5*la*n+3*N
98 +pow(x: la,y: 4)*N+6*la*la*N)
99 + pow(x: b2,y: 3)*(15+pow(x: la,y: 6)+15*pow(x: la,y: 4)+45*la*la)
100 + (pow(x: b3,y: 3)+3*b2*b2*b3+3*b3*b3*b2)
101 *(pow(x: la,y: 5)*n+14*pow(x: la,y: 3)*n+33*la*n+15*N
102 +15*pow(x: la,y: 4)*N+45*la*la*N+pow(x: la,y: 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(x: la,y: 4)+6*la*la); // ok
105 v1 = -2*b2*la - 2*b3*(n+la*N); // ok
106 v2 = -4*b2*b2*(3*la+pow(x: la,y: 3))
107 - (4*b3*b3+8*b2*b3)*(la*la*n+2*n+pow(x: la,y: 3)*N+3*la*N)
108 - 4*b1*b2*la - 4*b3*b1*(n+la*N); // ok
109 v3 = -12*b3*b1*(b3+2*b2)*(la*la*n+2*n+pow(x: la,y: 3)*N+3*la*N)
110 - 6*pow(x: b2,y: 3)*la*(15+pow(x: la,y: 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(x: la,y: 4)*n+pow(x: la,y: 5)*N
113 +10*pow(x: la,y: 3)*N)
114 - 6*b1*b1*b2*la - 6*b3*b1*b1*(n+la*N)
115 - 12*b2*b2*b1*(3*la+std::pow(x: la,y: 3)); // ok
116 z1 = b1 + b2*(3+la*la) + b3*(la*n+3*N+la*la*N); // ok
117 z2 = b1*b1 + b2*b2*(15+pow(x: la,y: 4)+18*la*la)
118 + (b3*b3+2*b2*b3)*(pow(x: la,y: 3)*n+17*la*n+15*N
119 +pow(x: la,y: 4)*N+18*la*la*N)
120 + 2*b1*b2*(3+la*la) + 2*b3*b1*(la*n+3*N+la*la*N); // ok
121 x1 = -6*b2*la - 2*b3*(4*n+3*la*N); // ok
122 b1_ = b1; b2_ = b2; b3_ = b3; la_ = la;
123 m1_ = m1; m2_ = m2; m3_ = m3;
124 v1_ = v1; v2_ = v2; v3_ = v3; z1_ = z1; z2_ = z2; x1_ = x1;
125 } else {
126 // these assignments are never used ?
127 // b1 = b1_; b2 = b2_; b3 = b3_; la = la_;
128 // m1 = m1_; m2 = m2_; m3 = m3_;
129 // v1 = v1_; v2 = v2_; v3 = v3_; z1 = z1_; z2 = z2_; x1 = x1_;
130 constants_match = true;
131 }
132
133 // compute the first four moments
134 if (!init_ || !constants_match || b0 != b0_ || h1 != h1_ || T != T_) {
135 // these assignments are never used ?
136 //b1 = b1_; b2 = b2_; b3 = b3_; la = la_;
137 m1 = m1_; m2 = m2_; m3 = m3_;
138 v1 = v1_; v2 = v2_; /*v3 = v3_;*/ z1 = z1_; /*z2 = z2_;*/ x1 = x1_;
139
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;
148 }
149
150 for (i = 0; i < T; ++i) {
151 Real m1i = m1ai[i];
152 Real m2i = m2ai[i];
153 Real m3i = m3ai[i];
154
155 Real m1im2i = m1i-m2i, m1im3i = m1i-m3i, m2im3i = m2i-m3i;
156 Real Eh = b0*(1-m1i)/(1-m1) + m1i*h1; // ko
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)
160 + m2i*h1*h1; // ko
161 Real Eh3 = pow(x: b0,y: 3)*(
162 (1-m3i)/(1-m3)
163 + 3*m2*((1-m3i)/(1-m3)-m2im3i/(m2-m3))/(1-m2)
164 + 3*m1*((1-m3i)/(1-m3)-m1im3i/(m1-m3))/(1-m1)
165 + 6*m1*m2*(
166 ((1-m3i)/(1-m3)-m2im3i/(m2-m3))/(1-m2)
167 + (m2im3i/(m2-m3)-m1im3i/(m1-m3))/(m1-m2)
168 )/(1-m1))
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)
172 + m3i*h1*h1*h1; // ko
173 Real Eh3_2 = .375*std::pow(x: Eh,y: -0.5)*Eh2+.625*std::pow(x: Eh,y: 1.5);
174 Real Eh5_2 = 1.875*std::pow(x: Eh,y: 0.5)*Eh2-.875*std::pow(x: Eh,y: 2.5);
175 sEh += Eh;
176 sEh2 += Eh2;
177 sEh3 += Eh3;
178 for (j = 0; j < T-i-1; ++j) {
179 Real Ehh = b0*Eh*(1-m1ai[j+1])/(1-m1)+ Eh2*m1ai[j+1]; // ko
180 Real Ehh2 = b0*b0*Eh*((1+m1)*(1-m2ai[j+1])/(1-m2)
181 - 2*m1*(m1ai[j+1]
182 -m2ai[j+1])/(m1-m2))/(1-m1)
183 + 2*b0*m1*Eh2*(m1ai[j+1]-m2ai[j+1])/(m1-m2)
184 + m2ai[j+1]*Eh3; // ko
185 Real Eh2h = b0*Eh2*(1-m1ai[j+1])/(1-m1)
186 + m1ai[j+1]*Eh3; // ok
187 Real Eh1_2eh = v1*m1ai[j]*Eh3_2; // ko
188 Real Eh1_2eh2 = 2*b0*v1*(m1ai[j+1]
189 -m2ai[j+1])*Eh3_2/(m1-m2)
190 + v2*m2ai[j]*Eh5_2; // ko
191 Real Ehij = b0*(1-m1ai[i+j+1])/(1-m1)
192 + m1ai[i+j+1]*h1; // ko
193 Real Ehh3_2 = 0.375*Ehh2/std::sqrt(x: Ehij)
194 + 0.75*std::sqrt(x: Ehij)*Ehh
195 - 0.125*std::pow(x: Ehij,y: 1.5)*Eh; // ko
196 Real Eh3_2eh = v1*m1ai[j]*Eh5_2; // ko
197 Real Eh3_2e3h = x1*m1ai[j]*Eh5_2; // ok
198 Real Eh1_2eh3_2 = 0.375*Eh1_2eh2/std::sqrt(x: Ehij)
199 + 0.75*std::sqrt(x: Ehij)*Eh1_2eh; // ko
200 sEhh += Ehh;
201 sEh1_2eh += Eh1_2eh;
202 sEhh2 += Ehh2;
203 sEh2h += Eh2h;
204 sEh1_2eh2 += Eh1_2eh2;
205 sEh3_2eh += Eh3_2eh;
206 sEhe2h += b0*Eh*(1-m1ai[j+1])/(1-m1)
207 + z1*m1ai[j]*Eh2; // ko
208 sEh3_2e3h += Eh3_2e3h; // ok
209 for (k = 0; k < T-i-j-2; ++k) {
210 Real Ehhh = b0*Ehh*(1-m1ai[k+1])/(1-m1)
211 + m1ai[k+1]*Ehh2; //ko
212 Real Eh1_2ehh = b0*Eh1_2eh*(1-m1ai[k+1])/(1-m1)
213 + m1ai[k+1]*Eh1_2eh2; // ko
214 sEhhh += Ehhh;
215 sEh1_2ehh += Eh1_2ehh;
216 sEhh1_2eh += v1*m1ai[k]*Ehh3_2; // ko
217 sEh1_2eh1_2eh += v1*m1ai[k]*Eh1_2eh3_2; // ko
218 }
219 }
220 }
221
222 ex = T*r - 0.5*sEh;
223 SD1 = 2*sEhh + sEh2;
224 SD2 = sEh;
225 SD3 = sEh1_2eh;
226 ex2 = T*T*r*r - T*r*sEh + 0.25*SD1 + SD2 - SD3;
227 ST1 = 6*sEhhh + (3*sEhh2 + (3*sEh2h + sEh3));
228 ST2 = 3*sEh1_2eh;
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(x: T*r,y: 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(x: T*r,y: 4) - 2*pow(x: T*r,y: 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);
240
241 // compute variance, skewness, kurtosis
242 sigma = ex2 - ex*ex;
243 // 3rd central moment mu3
244 k3 = ex3 - 3*sigma*ex - ex*ex*ex;
245 // 4th central moment mu4
246 k4 = ex4 + 6*ex*ex*ex2 - 3*ex*ex*ex*ex - 4*ex*ex3;
247 k3 /= std::pow(x: sigma,y: 1.5); // 3rd standardized moment, ie skewness
248 k4 /= pow(x: sigma,y: 2); // 4th standardized moment, ie kurtosis
249 ex_ = ex; sigma_ = sigma;
250 k3_ = k3; k4_ = k4; r_ = r; T_ = T; b0_ = b0; h1_ = h1;
251 } else {
252 ex = ex_; sigma = sigma_;
253 k3 = k3_; k4 = k4_; r = r_; T = T_; /*b0 = b0_; h1 = h1_;*/ // never used ?
254 }
255
256 // compute call option price
257 stdev = std::sqrt(x: sigma);
258 del = (ex - r*T + sigma/2)/stdev;
259 d = (std::log(x: s/x) + (r*T+sigma/2))/stdev;
260 d_ = d+del;
261 C = s*std::exp(x: del*stdev)*CumulativeNormalDistribution()(d_)
262 - x*std::exp(x: -r*T)*CumulativeNormalDistribution()(d_-stdev);
263 A3 = s*std::exp(x: del*stdev)*stdev*((2*stdev-d_)
264 *std::exp(x: -d_*d_/2)/std::sqrt(x: 2*M_PI)
265 +sigma*CumulativeNormalDistribution()(d_))/6;
266 A4 = s*std::exp(x: del*stdev)*stdev*(
267 (d_*d_-1-3*stdev*(d_-stdev))*exp(x: -d_*d_/2)/std::sqrt(x: 2*M_PI)
268 -sigma*stdev*CumulativeNormalDistribution()(d_))/24;
269 Capp = C + k3*A3 + (k4-3)*A4;
270 init_ = true;
271
272 switch (payoff->optionType()) {
273 case Option::Call:
274 results_.value = Capp;
275 break;
276 case Option::Put:
277 results_.value = Capp+strikePrice*riskFreeDiscount/dividendDiscount
278 -spotPrice;
279 break;
280 default:
281 QL_FAIL("unknown option type");
282 }
283 }
284}
285

source code of quantlib/ql/pricingengines/vanilla/analyticgjrgarchengine.cpp