QuantLib: a free/open-source library for quantitative finance
Fully annotated sources - version 1.32
Loading...
Searching...
No Matches
garch.cpp
1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2006 Joseph Wang
5 Copyright (C) 2012 Liquidnet Holdings, Inc.
6
7 This file is part of QuantLib, a free-software/open-source library
8 for financial quantitative analysts and developers - http://quantlib.org/
9
10 QuantLib is free software: you can redistribute it and/or modify it
11 under the terms of the QuantLib license. You should have received a
12 copy of the license along with this program; if not, please email
13 <quantlib-dev@lists.sf.net>. The license is also available online at
14 <http://quantlib.org/license.shtml>.
15
16 This program is distributed in the hope that it will be useful, but WITHOUT
17 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 FOR A PARTICULAR PURPOSE. See the license for more details.
19*/
20
21#include <ql/math/autocovariance.hpp>
22#include <ql/math/optimization/leastsquare.hpp>
23#include <ql/math/optimization/simplex.hpp>
24#include <ql/models/volatility/garch.hpp>
25#include <utility>
26
27namespace QuantLib {
28
29 namespace {
30
31 const Real tol_level = 1.0e-8;
32
33 class Garch11Constraint : public Constraint {
34 private:
35 class Impl final : public Constraint::Impl {
36 Real gammaLower_, gammaUpper_;
37 public:
38 Impl (Real gammaLower, Real gammaUpper)
39 : gammaLower_(gammaLower), gammaUpper_(gammaUpper) {}
40 bool test(const Array& x) const override {
41 QL_REQUIRE(x.size() >= 3, "size of parameters vector < 3");
42 return x[0] > 0 && x[1] >= 0 && x[2] >= 0
43 && x[1] + x[2] < gammaUpper_
44 && x[1] + x[2] >= gammaLower_;
45 }
46 };
47 public:
48 Garch11Constraint(Real gammaLower, Real gammaUpper)
49 : Constraint(ext::shared_ptr<Constraint::Impl>(
50 new Garch11Constraint::Impl(gammaLower, gammaUpper))) {}
51 };
52
53
54 class Garch11CostFunction : public CostFunction {
55 public:
56 explicit Garch11CostFunction (const std::vector<Volatility> &);
57 Real value(const Array& x) const override;
58 Array values(const Array& x) const override;
59 void gradient(Array& grad, const Array& x) const override;
60 Real valueAndGradient(Array& grad, const Array& x) const override;
61
62 private:
63 const std::vector<Volatility> &r2_;
64 };
65
66 Garch11CostFunction::Garch11CostFunction(
67 const std::vector<Volatility> &r2)
68 : r2_(r2) {}
69
70 Real Garch11CostFunction::value(const Array& x) const {
71 Real retval(0.0);
72 Real sigma2 = 0;
73 Real u2 = 0;
74 for (auto r2 : r2_) {
75 sigma2 = x[0] + x[1] * u2 + x[2] * sigma2;
76 u2 = r2;
77 retval += std::log(sigma2) + u2 / sigma2;
78 }
79 return retval / (2.0*r2_.size());
80 }
81
82 Array Garch11CostFunction::values(const Array& x) const {
83 Array retval (r2_.size());
84 Real sigma2 = 0;
85 Real u2 = 0;
86 Size i = 0;
87 for (auto r2 : r2_) {
88 sigma2 = x[0] + x[1] * u2 + x[2] * sigma2;
89 u2 = r2;
90 retval[i++] = (std::log(sigma2) + u2 / sigma2)/(2.0*r2_.size());
91 }
92 return retval;
93 }
94
95 void Garch11CostFunction::gradient(Array& grad, const Array& x) const {
96 std::fill (grad.begin(), grad.end(), 0.0);
97 Real sigma2 = 0;
98 Real u2 = 0;
99 Real sigma2prev = sigma2;
100 Real u2prev = u2;
101 Real norm = 2.0 * r2_.size();
102 for (auto r2 : r2_) {
103 sigma2 = x[0] + x[1] * u2 + x[2] * sigma2;
104 u2 = r2;
105 Real w = (sigma2 - u2) / (sigma2*sigma2);
106 grad[0] += w;
107 grad[1] += u2prev * w;
108 grad[2] += sigma2prev * w;
109 u2prev = u2;
110 sigma2prev = sigma2;
111 }
112 std::transform(grad.begin(), grad.end(), grad.begin(),
113 [=](Real x) -> Real { return x / norm; });
114 }
115
116 Real Garch11CostFunction::valueAndGradient(Array& grad,
117 const Array& x) const {
118 std::fill (grad.begin(), grad.end(), 0.0);
119 Real retval(0.0);
120 Real sigma2 = 0;
121 Real u2 = 0;
122 Real sigma2prev = sigma2;
123 Real u2prev = u2;
124 Real norm = 2.0 * r2_.size();
125 for (auto r2 : r2_) {
126 sigma2 = x[0] + x[1] * u2 + x[2] * sigma2;
127 u2 = r2;
128 retval += std::log(sigma2) + u2 / sigma2;
129 Real w = (sigma2 - u2) / (sigma2*sigma2);
130 grad[0] += w;
131 grad[1] += u2prev * w;
132 grad[2] += sigma2prev * w;
133 u2prev = u2;
134 sigma2prev = sigma2;
135 }
136 std::transform(grad.begin(), grad.end(), grad.begin(),
137 [=](Real x) -> Real { return x / norm; });
138 return retval / norm;
139 }
140
141
142 class FitAcfProblem : public LeastSquareProblem {
143 public:
144 FitAcfProblem(Real A2, Array acf, std::vector<std::size_t> idx);
145 Size size() override;
146 void targetAndValue(const Array& x, Array& target, Array& fct2fit) override;
147 void targetValueAndGradient(const Array& x,
148 Matrix& grad_fct2fit,
149 Array& target,
150 Array& fct2fit) override;
151
152 private:
153 Real A2_;
154 Array acf_;
155 std::vector<std::size_t> idx_;
156 };
157
158 FitAcfProblem::FitAcfProblem(Real A2, Array acf, std::vector<std::size_t> idx)
159 : A2_(A2), acf_(std::move(acf)), idx_(std::move(idx)) {}
160
161 Size FitAcfProblem::size() { return idx_.size(); }
162
163 void FitAcfProblem::targetAndValue(const Array& x, Array& target,
164 Array& fct2fit) {
165 Real A4 = acf_[0] + A2_*A2_;
166 Real gamma = x[0];
167 Real beta = x[1];
168 target[0] = A2_*A2_/A4;
169 fct2fit[0] =
170 (1 - 3*gamma*gamma - 2*beta*beta + 4*beta*gamma)
171 / (3*(1 - gamma*gamma));
172 target[1] = acf_[1] / A4;
173 fct2fit[1] = gamma * (1 - fct2fit[0]) - beta;
174 for (std::size_t i = 2; i < idx_.size(); ++i) {
175 target[i] = acf_[idx_[i]] / A4;
176 fct2fit[i] = std::pow(gamma, (int)idx_[i]-1)* fct2fit[1];
177 }
178 }
179
180 void FitAcfProblem::targetValueAndGradient(const Array& x,
181 Matrix& grad_fct2fit,
182 Array& target,
183 Array& fct2fit) {
184 Real A4 = acf_[0] + A2_*A2_;
185 Real gamma = x[0];
186 Real beta = x[1];
187 target[0] = A2_*A2_/A4;
188 Real w1 = (1 - 3*gamma*gamma - 2*beta*beta + 4*beta*gamma);
189 Real w2 = (1 - gamma*gamma);
190 fct2fit[0] = w1 / (3*w2);
191 grad_fct2fit[0][0] = (2.0/3.0) * ((2*beta-3*gamma)*w2 + 2*w1*gamma) / (w2*w2);
192 grad_fct2fit[0][1] = (4.0/3.0) * (gamma - beta) / w2;
193 target[1] = acf_[1] / A4;
194 fct2fit[1] = gamma * (1 - fct2fit[0]) - beta;
195 grad_fct2fit[1][0] = (1 - fct2fit[0]) - gamma * grad_fct2fit[0][0];
196 grad_fct2fit[1][1] = -gamma * grad_fct2fit[0][1] - 1;
197 for (std::size_t i = 2; i < idx_.size(); ++i) {
198 target[i] = acf_[idx_[i]] / A4;
199 w1 = std::pow(gamma, (int)idx_[i]-1);
200 fct2fit[i] = w1 * fct2fit[1];
201 grad_fct2fit[i][0] = (idx_[i]-1) * (w1/gamma)*fct2fit[1] + w1*grad_fct2fit[1][0];
202 grad_fct2fit[i][1] = w1 * grad_fct2fit[1][1];
203 }
204 }
205
206
207 class FitAcfConstraint : public Constraint {
208 private:
209 class Impl final : public Constraint::Impl {
210 Real gammaLower_, gammaUpper_;
211 public:
212 Impl(Real gammaLower, Real gammaUpper)
213 : gammaLower_(gammaLower), gammaUpper_(gammaUpper) {}
214 bool test(const Array& x) const override {
215 QL_REQUIRE(x.size() >= 2, "size of parameters vector < 2");
216 return x[0] >= gammaLower_ && x[0] < gammaUpper_
217 && x[1] >= 0 && x[1] <= x[0];
218 }
219 };
220 public:
221 FitAcfConstraint(Real gammaLower, Real gammaUpper)
222 : Constraint(ext::shared_ptr<Constraint::Impl>(
223 new FitAcfConstraint::Impl(gammaLower, gammaUpper))) {}
224 };
225
226
227 // Initial guess based on fitting ACF - initial guess for
228 // fitting acf is a moment matching estimates for mean(r2),
229 // acf(0), and acf(1).
230 Real initialGuess1(const Array &acf, Real mean_r2,
231 Real &alpha, Real &beta, Real &omega) {
232 Real A21 = acf[1];
233 Real A4 = acf[0] + mean_r2*mean_r2;
234
235 Real A = mean_r2*mean_r2/A4; // 1/sigma^2
236 Real B = A21 / A4; // rho(1)
237
238 Real gammaLower = A <= 1./3. - tol_level ? std::sqrt((1 - 3*A)/(3 - 3*A)) + tol_level : Real(tol_level);
239 Garch11Constraint constraints(gammaLower, 1.0 - tol_level);
240
241 Real gamma = gammaLower + (1 - gammaLower) * 0.5;
242 beta = std::min(gamma, std::max(gamma * (1 - A) - B, 0.0));
243 alpha = gamma - beta;
244 omega = mean_r2 * (1 - gamma);
245
246 if (std::fabs(A-0.5) < QL_EPSILON) {
247 gamma = std::max(gammaLower, -(1+4*B*B)/(4*B));
248 beta = std::min(gamma, std::max(gamma * (1 - A) - B, 0.0));
249 alpha = gamma - beta;
250 omega = mean_r2 * (1 - gamma);
251 } else {
252 if (A > 1.0 - QL_EPSILON) {
253 gamma = std::max(gammaLower, -(1+B*B)/(2*B));
254 beta = std::min(gamma, std::max(gamma * (1 - A) - B, 0.0));
255 alpha = gamma - beta;
256 omega = mean_r2 * (1 - gamma);
257 } else {
258 Real D = (3*A-1)*(2*B*B+(1-A)*(2*A-1));
259 if (D >= 0) {
260 Real d = std::sqrt(D);
261 Real b = (B - d)/(2*A-1);
262 Real g = 0;
263 if (b >= tol_level && b <= 1.0 - tol_level) {
264 g = (b + B) / (1 - A);
265 }
266 if (g < gammaLower) {
267 b = (B + d)/(2*A-1);
268 if (b >= tol_level && b <= 1.0 - tol_level) {
269 g = (b + B) / (1 - A);
270 }
271 }
272 if (g >= gammaLower) {
273 gamma = g;
274 beta = std::min(gamma, std::max(gamma * (1 - A) - B, 0.0));
275 alpha = gamma - beta;
276 omega = mean_r2 * (1 - gamma);
277 }
278 }
279 }
280 }
281
282 std::vector<std::size_t> idx;
283 std::size_t nCov = acf.size() - 1;
284 for (std::size_t i = 0; i <= nCov; ++i) {
285 if (i < 2 || (acf[i] > 0 && acf[i-1] > 0 && acf[i-1] > acf[i])) {
286 idx.push_back(i);
287 }
288 }
289
290 Array x(2);
291 x[0] = gamma;
292 x[1] = beta;
293
294 try {
295 FitAcfConstraint c(gammaLower, 1.0 - tol_level);
296 NonLinearLeastSquare nnls(c);
297 nnls.setInitialValue(x);
298 FitAcfProblem pr(mean_r2, acf, idx);
299 x = nnls.perform(pr);
300 Array guess(3);
301 guess[0] = mean_r2 * (1 - x[0]);
302 guess[1] = x[0] - x[1];
303 guess[2] = x[1];
304 if (constraints.test(guess)) {
305 omega = guess[0];
306 alpha = guess[1];
307 beta = guess[2];
308 }
309 } catch (const std::exception &) {
310 // failed -- returning initial values
311 }
312 return gammaLower;
313 }
314
315 // Initial guess based on fitting ACF - initial guess for
316 // fitting acf is an estimate of gamma = alpfa+beta based on
317 // the property: acf(i+1) = gamma*acf(i) for i > 1.
318 Real initialGuess2 (const Array &acf, Real mean_r2,
319 Real &alpha, Real &beta, Real &omega) {
320 Real A21 = acf[1];
321 Real A4 = acf[0] + mean_r2*mean_r2;
322 Real A = mean_r2*mean_r2/A4; // 1/sigma^2
323 Real B = A21 / A4; // rho(1)
324 Real gammaLower = A <= 1./3. - tol_level ? std::sqrt((1 - 3*A)/(3 - 3*A)) + tol_level : Real(tol_level);
325 Garch11Constraint constraints(gammaLower, 1.0 - tol_level);
326
327 // ACF
328 Real gamma = 0;
329 std::size_t nn = 0;
330 std::vector<std::size_t> idx;
331 std::size_t nCov = acf.size() - 1;
332 for (std::size_t i = 0; i <= nCov; ++i) {
333 if (i < 2) idx.push_back(i);
334 if (i > 1 && acf[i] > 0 && acf[i-1] > 0 && acf[i-1] > acf[i]) {
335 gamma += acf[i]/acf[i-1];
336 nn++;
337 idx.push_back(i);
338 }
339 }
340 if (nn > 0)
341 gamma /= nn;
342 if (gamma < gammaLower) gamma = gammaLower;
343 beta = std::min(gamma, std::max(gamma * (1 - A) - B, 0.0));
344 omega = mean_r2 * (1 - gamma);
345
346 Array x(2);
347 x[0] = gamma;
348 x[1] = beta;
349
350 try {
351 FitAcfConstraint c(gammaLower, 1 - tol_level);
352 NonLinearLeastSquare nnls(c);
353 nnls.setInitialValue(x);
354 FitAcfProblem pr(mean_r2, acf, idx);
355 x = nnls.perform(pr);
356 Array guess(3);
357 guess[0] = mean_r2 * (1 - x[0]);
358 guess[1] = x[0] - x[1];
359 guess[2] = x[1];
360 if (constraints.test(guess)) {
361 omega = guess[0];
362 alpha = guess[1];
363 beta = guess[2];
364 }
365 } catch (const std::exception &) {
366 // failed -- returning initial values
367 }
368 return gammaLower;
369 }
370
371 }
372
373 Garch11::time_series
374 Garch11::calculate(const time_series& quoteSeries,
375 Real alpha, Real beta, Real omega) {
376 time_series retval;
377 auto cur = quoteSeries.cbegin();
378 Real u = cur->second;
379 Real sigma2 = u*u;
380 while (++cur != quoteSeries.end()) {
381 sigma2 = omega + alpha * u * u + beta * sigma2;
382 retval[cur->first] = std::sqrt(sigma2);
383 u = cur->second;
384 }
385 sigma2 = omega + alpha * u * u + beta * sigma2;
386 --cur;
387 auto prev = cur;
388 retval[cur->first + (cur->first - (--prev)->first) ] = std::sqrt(sigma2);
389 return retval;
390 }
391
392
393 ext::shared_ptr<Problem> Garch11::calibrate_r2(
394 Mode mode, const std::vector<Volatility> &r2, Real mean_r2,
395 Real &alpha, Real &beta, Real &omega) {
396 EndCriteria endCriteria(10000, 500, tol_level, tol_level, tol_level);
397 Simplex method(0.001);
398 return calibrate_r2(mode, r2, mean_r2, method, endCriteria,
399 alpha, beta, omega);
400 }
401
402 ext::shared_ptr<Problem> Garch11::calibrate_r2(
403 Mode mode, const std::vector<Volatility> &r2, Real mean_r2,
404 OptimizationMethod &method, const EndCriteria &endCriteria,
405 Real &alpha, Real &beta, Real &omega) {
406 Real dataSize = Real(r2.size());
407 alpha = 0.0;
408 beta = 0.0;
409 omega = 0.0;
410 QL_REQUIRE (dataSize >= 4,
411 "Data series is too short to fit GARCH model");
412 QL_REQUIRE (mean_r2 > 0, "Data series is constant");
413 omega = mean_r2 * dataSize / (dataSize - 1);
414
415 // ACF
416 Size maxLag = (Size)std::sqrt(dataSize);
417 Array acf(maxLag+1);
418 std::vector<Volatility> tmp(r2.size());
419 std::transform (r2.begin(), r2.end(), tmp.begin(),
420 [=](Real x) -> Real { return x - mean_r2; });
421 autocovariances (tmp.begin(), tmp.end(), acf.begin(), maxLag);
422 QL_REQUIRE (acf[0] > 0, "Data series is constant");
423
424 Garch11CostFunction cost (r2);
425
426 // two initial guesses based on fitting ACF
427 Real gammaLower = 0.0;
428 Array opt1(3);
429 Real fCost1 = QL_MAX_REAL;
430 if (mode != GammaGuess) {
431 gammaLower = initialGuess1(acf, mean_r2, opt1[1], opt1[2], opt1[0]);
432 fCost1 = cost.value(opt1);
433 }
434
435 Array opt2(3);
436 Real fCost2 = QL_MAX_REAL;
437 if (mode != MomentMatchingGuess) {
438 gammaLower = initialGuess2(acf, mean_r2, opt2[1], opt2[2], opt2[0]);
439 fCost2 = cost.value(opt2);
440 }
441
442 Garch11Constraint constraints(gammaLower, 1.0 - tol_level);
443
444 ext::shared_ptr<Problem> ret;
445 if (mode != DoubleOptimization) {
446 try {
447 ret = calibrate_r2(r2, method, constraints, endCriteria,
448 fCost1 <= fCost2 ? opt1 : opt2,
449 alpha, beta, omega);
450 } catch (const std::exception &) {
451 if (fCost1 <= fCost2) {
452 alpha = opt1[1];
453 beta = opt1[2];
454 omega = opt1[0];
455 } else {
456 alpha = opt2[1];
457 beta = opt2[2];
458 omega = opt2[0];
459 }
460 }
461 } else {
462 ext::shared_ptr<Problem> ret1, ret2;
463 try {
464 ret1 = calibrate_r2(r2, method, constraints, endCriteria,
465 opt1, alpha, beta, omega);
466 opt1[1] = alpha;
467 opt1[2] = beta;
468 opt1[0] = omega;
469 if (constraints.test(opt1))
470 fCost1 = std::min(fCost1, cost.value(opt1));
471 } catch (const std::exception &) {
472 fCost1 = QL_MAX_REAL;
473 }
474
475 try {
476 ret2 = calibrate_r2(r2, method, constraints, endCriteria,
477 opt2, alpha, beta, omega);
478 opt2[1] = alpha;
479 opt2[2] = beta;
480 opt2[0] = omega;
481 if (constraints.test(opt2))
482 fCost2 = std::min(fCost2, cost.value(opt2));
483 } catch (const std::exception &) {
484 fCost2 = QL_MAX_REAL;
485 }
486
487 if (fCost1 <= fCost2) {
488 alpha = opt1[1];
489 beta = opt1[2];
490 omega = opt1[0];
491 ret = ret1;
492 } else {
493 alpha = opt2[1];
494 beta = opt2[2];
495 omega = opt2[0];
496 ret = ret2;
497 }
498 }
499 return ret;
500 }
501
502 ext::shared_ptr<Problem> Garch11::calibrate_r2(
503 const std::vector<Volatility> &r2,
504 OptimizationMethod &method,
505 const EndCriteria &endCriteria,
506 const Array &initGuess, Real &alpha, Real &beta, Real &omega) {
507 Garch11Constraint constraints(0.0, 1.0 - tol_level);
508 return calibrate_r2(r2, method, constraints, endCriteria,
509 initGuess, alpha, beta, omega);
510 }
511
512 ext::shared_ptr<Problem> Garch11::calibrate_r2(
513 const std::vector<Volatility> &r2,
514 Real mean_r2,
515 OptimizationMethod &method,
516 const EndCriteria &endCriteria,
517 const Array &initGuess, Real &alpha, Real &beta, Real &omega) {
518 std::vector<Volatility> tmp(r2.size());
519 std::transform(r2.begin(), r2.end(), tmp.begin(),
520 [=](Real x) -> Real { return x - mean_r2; });
521 return calibrate_r2(tmp, method, endCriteria, initGuess,
522 alpha, beta, omega);
523 }
524
525 ext::shared_ptr<Problem> Garch11::calibrate_r2(
526 const std::vector<Volatility> &r2,
527 OptimizationMethod &method,
528 Constraint &constraints,
529 const EndCriteria &endCriteria,
530 const Array &initGuess, Real &alpha, Real &beta, Real &omega) {
531 Garch11CostFunction cost(r2);
532 ext::shared_ptr<Problem> problem(
533 new Problem(cost, constraints, initGuess));
534 // TODO: check return value from minimize()
535 /* EndCriteria::Type ret = */
536 method.minimize(*problem, endCriteria);
537 const Array &optimum = problem->currentValue();
538 alpha = optimum[1];
539 beta = optimum[2];
540 omega = optimum[0];
541 return problem;
542 }
543
544 ext::shared_ptr<Problem> Garch11::calibrate_r2(
545 const std::vector<Volatility> &r2,
546 Real mean_r2,
547 OptimizationMethod &method,
548 Constraint &constraints,
549 const EndCriteria &endCriteria,
550 const Array &initGuess, Real &alpha, Real &beta, Real &omega) {
551 std::vector<Volatility> tmp(r2.size());
552 std::transform(r2.begin(), r2.end(), tmp.begin(),
553 [=](Real x) -> Real { return x - mean_r2; });
554 return calibrate_r2(tmp, method, constraints, endCriteria,
555 initGuess, alpha, beta, omega);
556 }
557
558}
559
QuantLib::Array
1-D array used in linear algebra.
Definition: array.hpp:52
QuantLib::Array::begin
const_iterator begin() const
Definition: array.hpp:503
QuantLib::Constraint
Base constraint class.
Definition: constraint.hpp:35
QuantLib::EndCriteria
Criteria to end optimization process:
Definition: endcriteria.hpp:40
QuantLib::OptimizationMethod
Abstract class for constrained optimization method.
Definition: method.hpp:36
QuantLib::OptimizationMethod::minimize
virtual EndCriteria::Type minimize(Problem &P, const EndCriteria &endCriteria)=0
minimize the optimization problem P
QuantLib::Problem
Constrained optimization problem.
Definition: problem.hpp:42
QuantLib::Simplex
Multi-dimensional simplex class.
Definition: simplex.hpp:58
QuantLib::TimeSeries
Container for historical data.
Definition: timeseries.hpp:51
QuantLib::TimeSeries::cbegin
const_iterator cbegin() const
Definition: timeseries.hpp:318
QuantLib::TimeSeries::end
const_iterator end() const
Definition: timeseries.hpp:159
QL_MAX_REAL
#define QL_MAX_REAL
Definition: qldefines.hpp:176
QL_EPSILON
#define QL_EPSILON
Definition: qldefines.hpp:178
QuantLib::Real
QL_REAL Real
real number
Definition: types.hpp:50
QuantLib::Size
std::size_t Size
size of a container
Definition: types.hpp:58
QuantLib::exponential_integrals_helper::g
Real g(Real x)
Definition: exponentialintegrals.cpp:59
QuantLib
Definition: any.hpp:35
QuantLib::autocovariances
void autocovariances(ForwardIterator begin, ForwardIterator end, OutputIterator out, std::size_t maxLag)
Unbiased auto-covariances.
Definition: autocovariance.hpp:99
std
STL namespace.
Generated by Doxygen 1.9.5