1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2012 Liquidnet Holdings, Inc.
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 "garch.hpp"
21#include "utilities.hpp"
22#include <ql/models/volatility/garch.hpp>
23#include <ql/time/calendars/target.hpp>
24#include <ql/math/optimization/levenbergmarquardt.hpp>
25#include <ql/math/randomnumbers/inversecumulativerng.hpp>
26#include <ql/math/randomnumbers/mt19937uniformrng.hpp>
27#include <ql/math/distributions/normaldistribution.hpp>
28
29using namespace QuantLib;
30using boost::unit_test_framework::test_suite;
31
32namespace garch_test {
33
34 class DummyOptimizationMethod : public OptimizationMethod {
35 public:
36 EndCriteria::Type minimize(Problem& P, const EndCriteria& endCriteria) override {
37 P.setFunctionValue(P.value(x: P.currentValue()));
38 return EndCriteria::None;
39 }
40 };
41
42 struct Results {
43 Real alpha;
44 Real beta;
45 Real omega;
46 Real logLikelihood;
47 };
48
49 typedef InverseCumulativeRng<MersenneTwisterUniformRng,
50 InverseCumulativeNormal>
51 GaussianGenerator;
52}
53
54#define CHECK(results, garch, member) \
55 if (std::fabs(results.member - garch.member()) > 1.0e-6) { \
56 BOOST_ERROR("Failed to reproduce expected " #member \
57 << "\n calculated: " << garch.member() \
58 << "\n expected: " << results.member); \
59 }
60
61void GARCHTest::testCalibration() {
62
63 BOOST_TEST_MESSAGE("Testing GARCH model calibration...");
64
65 using namespace garch_test;
66
67 Date start(7, July, 1962), d = start;
68 TimeSeries<Volatility> ts;
69 Garch11 garch(0.2, 0.3, 0.4);
70 GaussianGenerator rng(MersenneTwisterUniformRng(48));
71
72 Volatility r = 0.0, v = 0.0;
73 for (std::size_t i = 0; i < 50000; ++i, d += 1) {
74 v = garch.forecast(r, sigma2: v);
75 r = rng.next().value * std::sqrt(x: v);
76 ts[d] = r;
77 }
78
79 // Default calibration; works fine in most cases
80 Garch11 cgarch1(ts);
81
82 Results calibrated = { .alpha: 0.207592, .beta: 0.281979, .omega: 0.204647, .logLikelihood: -0.0217413 };
83
84 CHECK(calibrated, cgarch1, alpha);
85 CHECK(calibrated, cgarch1, beta);
86 CHECK(calibrated, cgarch1, omega);
87 CHECK(calibrated, cgarch1, logLikelihood);
88
89 // Type 1 initial guess - no further optimization
90 Garch11 cgarch2(ts, Garch11::MomentMatchingGuess);
91 DummyOptimizationMethod m;
92 cgarch2.calibrate(quoteSeries: ts, method&: m, endCriteria: EndCriteria (3, 2, 0.0, 0.0, 0.0));
93 Results expected1 = { .alpha: 0.265749, .beta: 0.156956, .omega: 0.230964, .logLikelihood: -0.0227179 };
94
95 CHECK(expected1, cgarch2, alpha);
96 CHECK(expected1, cgarch2, beta);
97 CHECK(expected1, cgarch2, omega);
98 CHECK(expected1, cgarch2, logLikelihood);
99
100 // Optimization from this initial guess
101 cgarch2.calibrate(quoteSeries: ts);
102
103 CHECK(calibrated, cgarch2, alpha);
104 CHECK(calibrated, cgarch2, beta);
105 CHECK(calibrated, cgarch2, omega);
106 CHECK(calibrated, cgarch2, logLikelihood);
107
108 // Type 2 initial guess - no further optimization
109 Garch11 cgarch3(ts, Garch11::GammaGuess);
110 cgarch3.calibrate(quoteSeries: ts, method&: m, endCriteria: EndCriteria (3, 2, 0.0, 0.0, 0.0));
111 Results expected2 = { .alpha: 0.269896, .beta: 0.211373, .omega: 0.207534, .logLikelihood: -0.022798 };
112
113 CHECK(expected2, cgarch3, alpha);
114 CHECK(expected2, cgarch3, beta);
115 CHECK(expected2, cgarch3, omega);
116 CHECK(expected2, cgarch3, logLikelihood);
117
118 // Optimization from this initial guess
119 cgarch3.calibrate(quoteSeries: ts);
120
121 CHECK(calibrated, cgarch3, alpha);
122 CHECK(calibrated, cgarch3, beta);
123 CHECK(calibrated, cgarch3, omega);
124 CHECK(calibrated, cgarch3, logLikelihood);
125
126 // Double optimization using type 1 and 2 initial guesses
127 Garch11 cgarch4(ts, Garch11::DoubleOptimization);
128 cgarch4.calibrate(quoteSeries: ts);
129
130 CHECK(calibrated, cgarch4, alpha);
131 CHECK(calibrated, cgarch4, beta);
132 CHECK(calibrated, cgarch4, omega);
133 CHECK(calibrated, cgarch4, logLikelihood);
134
135 // Alternative, gradient based optimization - usually gives worse
136 // results than simplex
137 LevenbergMarquardt lm;
138 cgarch4.calibrate(quoteSeries: ts, method&: lm, endCriteria: EndCriteria (100000, 500, 1e-8, 1e-8, 1e-8));
139 Results expected3 = { .alpha: 0.265196, .beta: 0.277364, .omega: 0.678812, .logLikelihood: -0.216313 };
140
141 CHECK(expected3, cgarch4, alpha);
142 CHECK(expected3, cgarch4, beta);
143 CHECK(expected3, cgarch4, omega);
144 CHECK(expected3, cgarch4, logLikelihood);
145}
146
147namespace garch_test {
148
149 static Real expected_calc[] = {
150 0.452769, 0.513323, 0.530141, 0.5350841, 0.536558,
151 0.536999, 0.537132, 0.537171, 0.537183, 0.537187
152 };
153
154 void check_ts(const std::pair<Date, Volatility> &x) {
155 if (x.first.serialNumber() < 22835 || x.first.serialNumber() > 22844) {
156 BOOST_ERROR("Failed to reproduce calculated GARCH time: "
157 << "\n calculated: " << x.first.serialNumber()
158 << "\n expected: [22835, 22844]");
159 }
160 Real error =
161 std::fabs(x: x.second - expected_calc[x.first.serialNumber()-22835]);
162 if (error > 1.0e-6) {
163 BOOST_ERROR("Failed to reproduce calculated GARCH value at "
164 << x.first.serialNumber() << ": "
165 << "\n calculated: " << x.second
166 << "\n expected: "
167 << expected_calc[x.first.serialNumber()-22835]);
168 }
169 }
170
171}
172
173void GARCHTest::testCalculation() {
174 BOOST_TEST_MESSAGE("Testing GARCH model calculation...");
175
176 using namespace garch_test;
177
178 Date d(7, July, 1962);
179 TimeSeries<Volatility> ts;
180 Garch11 garch(0.2, 0.3, 0.4);
181
182 Volatility r = 0.1;
183 for (std::size_t i = 0; i < 10; ++i, d += 1) {
184 ts[d] = r;
185 }
186
187 TimeSeries<Volatility> tsout = garch.calculate(quoteSeries: ts);
188 std::for_each(first: tsout.cbegin(), last: tsout.cend(), f: check_ts);
189}
190
191test_suite* GARCHTest::suite() {
192 auto* suite = BOOST_TEST_SUITE("GARCH model tests");
193 suite->add(QUANTLIB_TEST_CASE(&GARCHTest::testCalibration));
194 suite->add(QUANTLIB_TEST_CASE(&GARCHTest::testCalculation));
195 return suite;
196}
197

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