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 <ql/math/distributions/chisquaredistribution.hpp>
21#include <ql/math/distributions/normaldistribution.hpp>
22#include <ql/processes/eulerdiscretization.hpp>
23#include <ql/processes/gjrgarchprocess.hpp>
24#include <ql/quotes/simplequote.hpp>
25#include <utility>
26
27namespace QuantLib {
28
29 GJRGARCHProcess::GJRGARCHProcess(Handle<YieldTermStructure> riskFreeRate,
30 Handle<YieldTermStructure> dividendYield,
31 Handle<Quote> s0,
32 Real v0,
33 Real omega,
34 Real alpha,
35 Real beta,
36 Real gamma,
37 Real lambda,
38 Real daysPerYear,
39 Discretization d)
40 : StochasticProcess(ext::shared_ptr<discretization>(new EulerDiscretization)),
41 riskFreeRate_(std::move(riskFreeRate)), dividendYield_(std::move(dividendYield)),
42 s0_(std::move(s0)), v0_(v0), omega_(omega), alpha_(alpha), beta_(beta), gamma_(gamma),
43 lambda_(lambda), daysPerYear_(daysPerYear), discretization_(d) {
44 registerWith(h: riskFreeRate_);
45 registerWith(h: dividendYield_);
46 registerWith(h: s0_);
47 }
48
49 Size GJRGARCHProcess::size() const {
50 return 2;
51 }
52
53 Array GJRGARCHProcess::initialValues() const {
54 return { s0_->value(), daysPerYear_*v0_ };
55 }
56
57 Array GJRGARCHProcess::drift(Time t, const Array& x) const {
58 const Real N = CumulativeNormalDistribution()(lambda_);
59 const Real n = std::exp(x: -lambda_*lambda_/2.0)/std::sqrt(x: 2*M_PI);
60 const Real q2 = 1.0 + lambda_*lambda_;
61 const Real q3 = lambda_*n + N + lambda_*lambda_*N;
62 const Real vol = (x[1] > 0.0) ? std::sqrt(x: x[1])
63 : (discretization_ == Reflection) ? Real(-std::sqrt(x: -x[1]))
64 : 0.0;
65
66 return {
67 riskFreeRate_->forwardRate(t1: t, t2: t, comp: Continuous).rate()
68 - dividendYield_->forwardRate(t1: t, t2: t, comp: Continuous).rate()
69 - 0.5 * vol * vol,
70 daysPerYear_*daysPerYear_*omega_ + daysPerYear_*(beta_
71 + alpha_*q2 + gamma_*q3 - 1.0) *
72 ((discretization_==PartialTruncation) ? x[1] : vol*vol)
73 };
74 }
75
76 Matrix GJRGARCHProcess::diffusion(Time, const Array& x) const {
77 /* the correlation matrix is
78 | 1 rho |
79 | rho 1 |
80 whose square root (which is used here) is
81 | 1 0 |
82 | rho std::sqrt(1-rho^2) |
83 */
84 Matrix tmp(2,2);
85 const Real N = CumulativeNormalDistribution()(lambda_);
86 const Real n = std::exp(x: -lambda_*lambda_/2.0)/std::sqrt(x: 2*M_PI);
87 const Real sigma2 = 2.0 + 4.0*lambda_*lambda_;
88 const Real q3 = lambda_*n + N + lambda_*lambda_*N;
89 const Real Eml_e4 = lambda_*lambda_*lambda_*n + 5.0*lambda_*n
90 + 3.0*N + lambda_*lambda_*lambda_*lambda_*N
91 + 6.0*lambda_*lambda_*N;
92 const Real sigma3 = Eml_e4 - q3*q3;
93 const Real sigma12 = -2.0*lambda_;
94 const Real sigma13 = -2.0*n - 2*lambda_*N;
95 const Real sigma23 = 2.0*N + sigma12*sigma13;
96 const Real vol = (x[1] > 0.0) ? std::sqrt(x: x[1])
97 : (discretization_ == Reflection) ? Real(- std::sqrt(x: -x[1]))
98 : 1e-8; // set vol to (almost) zero but still
99 // expose some correlation information
100 const Real rho1 = std::sqrt(x: daysPerYear_)*(alpha_*sigma12
101 + gamma_*sigma13) * vol * vol;
102 const Real rho2 = vol*vol*std::sqrt(x: daysPerYear_)
103 *std::sqrt(x: alpha_*alpha_*(sigma2 - sigma12*sigma12)
104 + gamma_*gamma_*(sigma3 - sigma13*sigma13)
105 + 2.0*alpha_*gamma_*(sigma23 - sigma12*sigma13));
106
107 // tmp[0][0], tmp[0][1] are the coefficients of dW_1 and dW_2
108 // in asset return stochastic process
109 tmp[0][0] = vol; tmp[0][1] = 0.0;
110 tmp[1][0] = rho1; tmp[1][1] = rho2;
111 return tmp;
112 }
113
114 Array GJRGARCHProcess::apply(const Array& x0,
115 const Array& dx) const {
116 return { x0[0] * std::exp(x: dx[0]), x0[1] + dx[1] };
117 }
118
119 Array GJRGARCHProcess::evolve(Time t0, const Array& x0,
120 Time dt, const Array& dw) const {
121 Array retVal(2);
122 Real vol, mu, nu;
123
124 const Real sdt = std::sqrt(x: dt);
125 const Real N = CumulativeNormalDistribution()(lambda_);
126 const Real n = std::exp(x: -lambda_*lambda_/2.0)/std::sqrt(x: 2*M_PI);
127 const Real sigma2 = 2.0 + 4.0*lambda_*lambda_;
128 const Real q2 = 1.0 + lambda_*lambda_;
129 const Real q3 = lambda_*n + N + lambda_*lambda_*N;
130 const Real Eml_e4 = lambda_*lambda_*lambda_*n + 5.0*lambda_*n
131 + 3.0*N + lambda_*lambda_*lambda_*lambda_*N
132 + 6.0*lambda_*lambda_*N;
133 const Real sigma3 = Eml_e4 - q3*q3;
134 const Real sigma12 = -2.0*lambda_;
135 const Real sigma13 = -2.0*n - 2*lambda_*N;
136 const Real sigma23 = 2.0*N + sigma12*sigma13;
137 const Real rho1 = std::sqrt(x: daysPerYear_)*(alpha_*sigma12 + gamma_*sigma13);
138 const Real rho2 = std::sqrt(x: daysPerYear_)
139 *std::sqrt(x: alpha_*alpha_*(sigma2 - sigma12*sigma12)
140 + gamma_*gamma_*(sigma3 - sigma13*sigma13)
141 + 2.0*alpha_*gamma_*(sigma23 - sigma12*sigma13));
142
143 switch (discretization_) {
144 // For the definition of PartialTruncation, FullTruncation
145 // and Reflection see Lord, R., R. Koekkoek and D. van Dijk (2006),
146 // "A Comparison of biased simulation schemes for
147 // stochastic volatility models",
148 // Working Paper, Tinbergen Institute
149 case PartialTruncation:
150 vol = (x0[1] > 0.0) ? Real(std::sqrt(x: x0[1])) : 0.0;
151 mu = riskFreeRate_->forwardRate(t1: t0, t2: t0+dt, comp: Continuous).rate()
152 - dividendYield_->forwardRate(t1: t0, t2: t0+dt, comp: Continuous).rate()
153 - 0.5 * vol * vol;
154 nu = daysPerYear_*daysPerYear_*omega_
155 + daysPerYear_*(beta_ + alpha_*q2 + gamma_*q3 - 1.0) * x0[1];
156
157 retVal[0] = x0[0] * std::exp(x: mu*dt+vol*dw[0]*sdt);
158 retVal[1] = x0[1] + nu*dt + sdt*vol*vol*(rho1*dw[0] + rho2*dw[1]);
159 break;
160 case FullTruncation:
161 vol = (x0[1] > 0.0) ? Real(std::sqrt(x: x0[1])) : 0.0;
162 mu = riskFreeRate_->forwardRate(t1: t0, t2: t0+dt, comp: Continuous).rate()
163 - dividendYield_->forwardRate(t1: t0, t2: t0+dt, comp: Continuous).rate()
164 - 0.5 * vol * vol;
165 nu = daysPerYear_*daysPerYear_*omega_
166 + daysPerYear_*(beta_ + alpha_*q2 + gamma_*q3 - 1.0) * vol *vol;
167
168 retVal[0] = x0[0] * std::exp(x: mu*dt+vol*dw[0]*sdt);
169 retVal[1] = x0[1] + nu*dt + sdt*vol*vol*(rho1*dw[0] + rho2*dw[1]);
170 break;
171 case Reflection:
172 vol = std::sqrt(x: std::fabs(x: x0[1]));
173 mu = riskFreeRate_->forwardRate(t1: t0, t2: t0+dt, comp: Continuous).rate()
174 - dividendYield_->forwardRate(t1: t0, t2: t0+dt, comp: Continuous).rate()
175 - 0.5 * vol*vol;
176 nu = daysPerYear_*daysPerYear_*omega_
177 + daysPerYear_*(beta_ + alpha_*q2 + gamma_*q3 - 1.0) * vol * vol;
178
179 retVal[0] = x0[0]*std::exp(x: mu*dt+vol*dw[0]*sdt);
180 retVal[1] = vol*vol
181 +nu*dt + sdt*vol*vol*(rho1*dw[0] + rho2*dw[1]);
182 break;
183 default:
184 QL_FAIL("unknown discretization schema");
185 }
186
187 return retVal;
188 }
189
190 const Handle<Quote>& GJRGARCHProcess::s0() const {
191 return s0_;
192 }
193
194 const Handle<YieldTermStructure>& GJRGARCHProcess::dividendYield() const {
195 return dividendYield_;
196 }
197
198 const Handle<YieldTermStructure>& GJRGARCHProcess::riskFreeRate() const {
199 return riskFreeRate_;
200 }
201
202 Time GJRGARCHProcess::time(const Date& d) const {
203 return riskFreeRate_->dayCounter().yearFraction(
204 d1: riskFreeRate_->referenceDate(), d2: d);
205 }
206
207}
208

source code of quantlib/ql/processes/gjrgarchprocess.cpp