[go: up one dir, main page]
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#include "gjrgarchmodel.hpp"
21#include "utilities.hpp"
22#include <ql/processes/gjrgarchprocess.hpp>
23#include <ql/math/distributions/normaldistribution.hpp>
24#include <ql/models/equity/gjrgarchmodel.hpp>
25#include <ql/models/equity/hestonmodelhelper.hpp>
26#include <ql/pricingengines/vanilla/analyticgjrgarchengine.hpp>
27#include <ql/pricingengines/vanilla/mceuropeangjrgarchengine.hpp>
28#include <ql/pricingengines/blackformula.hpp>
29#include <ql/time/calendars/target.hpp>
30#include <ql/time/calendars/nullcalendar.hpp>
31#include <ql/time/daycounters/actual365fixed.hpp>
32#include <ql/time/daycounters/actual360.hpp>
33#include <ql/time/daycounters/actualactual.hpp>
34#include <ql/termstructures/yield/zerocurve.hpp>
35#include <ql/termstructures/yield/flatforward.hpp>
36#include <ql/math/optimization/simplex.hpp>
37#include <ql/time/period.hpp>
38#include <ql/quotes/simplequote.hpp>
39
40using namespace QuantLib;
41using namespace boost::unit_test_framework;
42
43void GJRGARCHModelTest::testEngines() {
44 BOOST_TEST_MESSAGE(
45 "Testing Monte Carlo GJR-GARCH engine against "
46 "analytic GJR-GARCH engine...");
47
48 DayCounter dayCounter = ActualActual(ActualActual::ISDA);
49
50 const Date today = Date::todaysDate();
51 Handle<YieldTermStructure> riskFreeTS(flatRate(today, forward: 0.05, dc: dayCounter));
52 Handle<YieldTermStructure> dividendTS(flatRate(today, forward: 0.0, dc: dayCounter));
53
54 const Real s0 = 50.0;
55 const Real omega = 2.0e-6;
56 const Real alpha = 0.024;
57 const Real beta = 0.93;
58 const Real gamma = 0.059;
59 const Real daysPerYear = 365.0; // number of trading days per year
60 const Size maturity[] = {90, 180};
61 const Real strike[] = {35,40,45,50,55,60};
62 const Real Lambda[] = {0.0,0.1,0.2};
63 Real analytic[3][2][6]; // correct values of analytic approximation
64 analytic[0][0][0] = 15.4315;
65 analytic[0][0][1] = 10.5552;
66 analytic[0][0][2] = 5.9625;
67 analytic[0][0][3] = 2.3282;
68 analytic[0][0][4] = 0.5408;
69 analytic[0][0][5] = 0.0835;
70 analytic[0][1][0] = 15.8969;
71 analytic[0][1][1] = 11.2173;
72 analytic[0][1][2] = 6.9112;
73 analytic[0][1][3] = 3.4788;
74 analytic[0][1][4] = 1.3769;
75 analytic[0][1][5] = 0.4357;
76 analytic[1][0][0] = 15.4556;
77 analytic[1][0][1] = 10.6929;
78 analytic[1][0][2] = 6.2381;
79 analytic[1][0][3] = 2.6831;
80 analytic[1][0][4] = 0.7822;
81 analytic[1][0][5] = 0.1738;
82 analytic[1][1][0] = 16.0587;
83 analytic[1][1][1] = 11.5338;
84 analytic[1][1][2] = 7.3170;
85 analytic[1][1][3] = 3.9074;
86 analytic[1][1][4] = 1.7279;
87 analytic[1][1][5] = 0.6568;
88 analytic[2][0][0] = 15.8000;
89 analytic[2][0][1] = 11.2734;
90 analytic[2][0][2] = 7.0376;
91 analytic[2][0][3] = 3.6767;
92 analytic[2][0][4] = 1.5871;
93 analytic[2][0][5] = 0.5934;
94 analytic[2][1][0] = 16.9286;
95 analytic[2][1][1] = 12.3170;
96 analytic[2][1][2] = 8.0405;
97 analytic[2][1][3] = 4.6348;
98 analytic[2][1][4] = 2.3429;
99 analytic[2][1][5] = 1.0590;
100 Real mcValues[3][2][6]; // correct values of Monte Carlo
101 mcValues[0][0][0] = 15.4332;
102 mcValues[0][0][1] = 10.5453;
103 mcValues[0][0][2] = 5.9351;
104 mcValues[0][0][3] = 2.3521;
105 mcValues[0][0][4] = 0.5597;
106 mcValues[0][0][5] = 0.0776;
107 mcValues[0][1][0] = 15.8910;
108 mcValues[0][1][1] = 11.1772;
109 mcValues[0][1][2] = 6.8827;
110 mcValues[0][1][3] = 3.5096;
111 mcValues[0][1][4] = 1.4196;
112 mcValues[0][1][5] = 0.4502;
113 mcValues[1][0][0] = 15.4580;
114 mcValues[1][0][1] = 10.6433;
115 mcValues[1][0][2] = 6.2019;
116 mcValues[1][0][3] = 2.7513;
117 mcValues[1][0][4] = 0.8374;
118 mcValues[1][0][5] = 0.1706;
119 mcValues[1][1][0] = 15.9884;
120 mcValues[1][1][1] = 11.4139;
121 mcValues[1][1][2] = 7.3103;
122 mcValues[1][1][3] = 4.0497;
123 mcValues[1][1][4] = 1.8862;
124 mcValues[1][1][5] = 0.7322;
125 mcValues[2][0][0] = 15.6619;
126 mcValues[2][0][1] = 11.1263;
127 mcValues[2][0][2] = 7.0968;
128 mcValues[2][0][3] = 3.9152;
129 mcValues[2][0][4] = 1.8133;
130 mcValues[2][0][5] = 0.7010;
131 mcValues[2][1][0] = 16.5195;
132 mcValues[2][1][1] = 12.3181;
133 mcValues[2][1][2] = 8.6085;
134 mcValues[2][1][3] = 5.5700;
135 mcValues[2][1][4] = 3.3103;
136 mcValues[2][1][5] = 1.8053;
137
138 for (Size k = 0; k < 3; ++k) {
139 Real lambda = Lambda[k];
140 Real m1 = beta+(alpha+gamma*CumulativeNormalDistribution()(lambda))
141 *(1.0+lambda*lambda)+gamma*lambda*std::exp(x: -lambda*lambda/2.0)
142 /std::sqrt(x: 2.0*M_PI);
143 Real v0 = omega/(1.0-m1);
144 Handle<Quote> q(ext::shared_ptr<Quote>(new SimpleQuote(s0)));
145 ext::shared_ptr<GJRGARCHProcess> process(new GJRGARCHProcess(
146 riskFreeTS, dividendTS, q, v0, omega, alpha, beta, gamma, lambda, daysPerYear));
147 ext::shared_ptr<PricingEngine> engine1 =
148 MakeMCEuropeanGJRGARCHEngine<PseudoRandom>(process)
149 .withStepsPerYear(steps: 20)
150 .withAbsoluteTolerance(tolerance: 0.02)
151 .withSeed(seed: 1234);
152
153 ext::shared_ptr<PricingEngine> engine2(
154 new AnalyticGJRGARCHEngine(ext::make_shared<GJRGARCHModel>(
155 args&: process)));
156 for (Size i = 0; i < 2; ++i) {
157 for (Size j = 0; j < 6; ++j) {
158 Real x = strike[j];
159
160 ext::shared_ptr<StrikedTypePayoff> payoff(
161 new PlainVanillaPayoff(Option::Call, x));
162 Date exDate = today + maturity[i];
163 ext::shared_ptr<Exercise> exercise(
164 new EuropeanExercise(exDate));
165
166 VanillaOption option(payoff, exercise);
167
168 option.setPricingEngine(engine1);
169 Real calculated = option.NPV();
170
171 option.setPricingEngine(engine2);
172 Real expected = option.NPV();
173 Real tolerance = 7.5e-2;
174
175 if (std::fabs(x: expected - analytic[k][i][j]) > 2.0*tolerance) {
176 BOOST_ERROR("failed to match results from engines"
177 << "\n correct value: "
178 << analytic[k][i][j]
179 << "\n Analytic Approx.: "
180 << expected
181 << " +/- " << tolerance);
182 }
183 if (std::fabs(x: calculated-mcValues[k][i][j]) > 2.0*tolerance) {
184 BOOST_ERROR("failed to match results from engines"
185 << "\n correct value: "
186 << mcValues[k][i][j]
187 << "\n Monte Carlo: " << calculated
188 << " +/- " << tolerance);
189 }
190 }
191 }
192 }
193}
194
195
196void GJRGARCHModelTest::testDAXCalibration() {
197 /* this example is taken from A. Sepp
198 Pricing European-Style Options under Jump Diffusion Processes
199 with Stochstic Volatility: Applications of Fourier Transform
200 http://math.ut.ee/~spartak/papers/stochjumpvols.pdf
201 */
202
203 BOOST_TEST_MESSAGE(
204 "Testing GJR-GARCH model calibration using DAX volatility data...");
205
206 Date settlementDate(5, July, 2002);
207 Settings::instance().evaluationDate() = settlementDate;
208
209 DayCounter dayCounter = Actual365Fixed();
210 Calendar calendar = TARGET();
211
212 Integer t[] = { 13, 41, 75, 165, 256, 345, 524, 703 };
213 Rate r[] = { 0.0357,0.0349,0.0341,0.0355,0.0359,0.0368,0.0386,0.0401 };
214
215 std::vector<Date> dates;
216 std::vector<Rate> rates;
217 dates.push_back(x: settlementDate);
218 rates.push_back(x: 0.0357);
219 Size i;
220 for (i = 0; i < 8; ++i) {
221 dates.push_back(x: settlementDate + t[i]);
222 rates.push_back(x: r[i]);
223 }
224 Handle<YieldTermStructure> riskFreeTS(
225 ext::shared_ptr<YieldTermStructure>(
226 new ZeroCurve(dates, rates, dayCounter)));
227
228 Handle<YieldTermStructure> dividendTS(
229 flatRate(today: settlementDate, forward: 0.0, dc: dayCounter));
230
231 Volatility v[] =
232 { 0.6625,0.4875,0.4204,0.3667,0.3431,0.3267,0.3121,0.3121,
233 0.6007,0.4543,0.3967,0.3511,0.3279,0.3154,0.2984,0.2921,
234 0.5084,0.4221,0.3718,0.3327,0.3155,0.3027,0.2919,0.2889,
235 0.4541,0.3869,0.3492,0.3149,0.2963,0.2926,0.2819,0.2800,
236 0.4060,0.3607,0.3330,0.2999,0.2887,0.2811,0.2751,0.2775,
237 0.3726,0.3396,0.3108,0.2781,0.2788,0.2722,0.2661,0.2686,
238 0.3550,0.3277,0.3012,0.2781,0.2781,0.2661,0.2661,0.2681,
239 0.3428,0.3209,0.2958,0.2740,0.2688,0.2627,0.2580,0.2620,
240 0.3302,0.3062,0.2799,0.2631,0.2573,0.2533,0.2504,0.2544,
241 0.3343,0.2959,0.2705,0.2540,0.2504,0.2464,0.2448,0.2462,
242 0.3460,0.2845,0.2624,0.2463,0.2425,0.2385,0.2373,0.2422,
243 0.3857,0.2860,0.2578,0.2399,0.2357,0.2327,0.2312,0.2351,
244 0.3976,0.2860,0.2607,0.2356,0.2297,0.2268,0.2241,0.2320 };
245
246 Handle<Quote> s0(ext::shared_ptr<Quote>(new SimpleQuote(4468.17)));
247 Real strike[] = { 3400,3600,3800,4000,4200,4400,
248 4500,4600,4800,5000,5200,5400,5600 };
249
250 std::vector<ext::shared_ptr<CalibrationHelper> > options;
251
252 const Real omega = 2.0e-6;
253 const Real alpha = 0.024;
254 const Real beta = 0.93;
255 const Real gamma = 0.059;
256 const Real lambda = 0.1;
257 const Real daysPerYear = 365.0; // number of trading days per year
258 const Real m1 = beta+(alpha+gamma*CumulativeNormalDistribution()(lambda))
259 *(1.0+lambda*lambda)+gamma*lambda*std::exp(x: -lambda*lambda/2.0)
260 /std::sqrt(x: 2.0*M_PI);
261 const Real v0 = omega/(1.0-m1);
262
263 ext::shared_ptr<GJRGARCHProcess> process(new GJRGARCHProcess(
264 riskFreeTS, dividendTS, s0, v0,
265 omega, alpha, beta, gamma, lambda, daysPerYear));
266 ext::shared_ptr<GJRGARCHModel> model(new GJRGARCHModel(process));
267
268 ext::shared_ptr<PricingEngine> engine(
269 new AnalyticGJRGARCHEngine(ext::shared_ptr<GJRGARCHModel>(model)));
270
271 for (Size s = 3; s < 10; ++s) {
272 for (Size m = 0; m < 3; ++m) {
273 Handle<Quote> vol(ext::shared_ptr<Quote>(
274 new SimpleQuote(v[s*8+m])));
275
276 Period maturity((int)((t[m]+3)/7.), Weeks); // round to weeks
277 ext::shared_ptr<BlackCalibrationHelper> option(
278 new HestonModelHelper(maturity, calendar,
279 s0->value(), strike[s], vol,
280 riskFreeTS, dividendTS,
281 BlackCalibrationHelper::ImpliedVolError));
282 option->setPricingEngine(engine);
283 options.push_back(x: option);
284 }
285 }
286
287 Simplex om(0.05);
288 model->calibrate(options, method&: om,
289 endCriteria: EndCriteria(400, 40, 1.0e-8, 1.0e-8, 1.0e-8));
290
291 Real sse = 0;
292 for (i = 0; i < options.size(); ++i) {
293 const Real diff = options[i]->calibrationError()*100.0;
294 sse += diff*diff;
295 }
296 Real maxExpected = 15;
297 if (sse > maxExpected) {
298 BOOST_FAIL("Failed to reproduce calibration error"
299 << "\n calculated: " << sse
300 << "\n expected: < " << maxExpected);
301 }
302}
303
304test_suite* GJRGARCHModelTest::suite(SpeedLevel speed) {
305 auto* suite = BOOST_TEST_SUITE("GJR-GARCH model tests");
306
307 if (speed <= Fast) {
308 suite->add(QUANTLIB_TEST_CASE(&GJRGARCHModelTest::testDAXCalibration));
309 }
310
311 if (speed == Slow) {
312 suite->add(QUANTLIB_TEST_CASE(&GJRGARCHModelTest::testEngines));
313 }
314
315 return suite;
316}
317

source code of quantlib/test-suite/gjrgarchmodel.cpp