1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2007, 2008 Klaus Spanderen
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 hybridhestonhullwhiteprocess.hpp
21 \brief hybrid equity (heston model)
22 with stochastic interest rates (hull white model)
23*/
24
25#include <ql/termstructures/yield/flatforward.hpp>
26#include <ql/processes/hybridhestonhullwhiteprocess.hpp>
27
28namespace QuantLib {
29
30 HybridHestonHullWhiteProcess::HybridHestonHullWhiteProcess(
31 const ext::shared_ptr<HestonProcess> & hestonProcess,
32 const ext::shared_ptr<HullWhiteForwardProcess> & hullWhiteProcess,
33 Real corrEquityShortRate,
34 HybridHestonHullWhiteProcess::Discretization discretization)
35 : hestonProcess_(hestonProcess),
36 hullWhiteProcess_(hullWhiteProcess),
37 hullWhiteModel_(new HullWhite(hestonProcess->riskFreeRate(),
38 hullWhiteProcess->a(),
39 hullWhiteProcess->sigma())),
40 corrEquityShortRate_(corrEquityShortRate),
41 discretization_(discretization),
42 maxRho_(std::sqrt(x: 1-hestonProcess->rho()*hestonProcess->rho())
43 - std::sqrt(QL_EPSILON) /* reserve for rounding errors */),
44
45 T_(hullWhiteProcess->getForwardMeasureTime()),
46 endDiscount_(hestonProcess->riskFreeRate()->discount(t: T_)) {
47
48 QL_REQUIRE( corrEquityShortRate*corrEquityShortRate
49 +hestonProcess->rho()*hestonProcess->rho() <= 1.0,
50 "correlation matrix is not positive definite");
51
52 QL_REQUIRE(hullWhiteProcess->sigma() > 0.0,
53 "positive vol of Hull White process is required");
54 }
55
56 Size HybridHestonHullWhiteProcess::size() const {
57 return 3;
58 }
59
60 Array HybridHestonHullWhiteProcess::initialValues() const {
61 return {
62 hestonProcess_->s0()->value(),
63 hestonProcess_->v0(),
64 hullWhiteProcess_->x0()
65 };
66 }
67
68 Array HybridHestonHullWhiteProcess::drift(Time t, const Array& x) const {
69 Array x0 = { x[0], x[1] };
70 Array y0 = hestonProcess_->drift(t, x: x0);
71
72 return {
73 y0[0],
74 y0[1],
75 hullWhiteProcess_->drift(t, x: x[2])
76 };
77 }
78
79 Array HybridHestonHullWhiteProcess::apply(const Array& x0, const Array& dx) const {
80 Array xt = { x0[0], x0[1] }, dxt = { dx[0], dx[1] };
81 Array yt = hestonProcess_->apply(x0: xt, dx: dxt);
82
83 return {
84 yt[0],
85 yt[1],
86 hullWhiteProcess_->apply(x0: x0[2], dx: dx[2])
87 };
88 }
89
90 Matrix HybridHestonHullWhiteProcess::diffusion(Time t, const Array& x) const {
91 Matrix retVal(3,3);
92
93 Array xt(2); xt[0] = x[0]; xt[1] = x[1];
94 Matrix m = hestonProcess_->diffusion(t, x: xt);
95 retVal[0][0] = m[0][0]; retVal[0][1] = 0.0; retVal[0][2] = 0.0;
96 retVal[1][0] = m[1][0]; retVal[1][1] = m[1][1]; retVal[1][2] = 0.0;
97
98 const Real sigma = hullWhiteProcess_->sigma();
99 retVal[2][0] = corrEquityShortRate_ * sigma;
100 retVal[2][1] = - retVal[2][0]*retVal[1][0] / retVal[1][1];
101 retVal[2][2] = std::sqrt( x: sigma*sigma - retVal[2][1]*retVal[2][1]
102 - retVal[2][0]*retVal[2][0] );
103
104 return retVal;
105 }
106
107 Array HybridHestonHullWhiteProcess::evolve(Time t0, const Array& x0,
108 Time dt, const Array& dw) const {
109
110 const Rate r = x0[2];
111 const Real a = hullWhiteProcess_->a();
112 const Real sigma = hullWhiteProcess_->sigma();
113 const Real rho = corrEquityShortRate_;
114 const Real xi = hestonProcess_->rho();
115 const Volatility eta = (x0[1] > 0.0) ? Real(std::sqrt(x: x0[1])) : 0.0;
116 const Time s = t0;
117 const Time t = t0 + dt;
118 const Time T = T_;
119 const Rate dy
120 = hestonProcess_->dividendYield()->forwardRate(t1: s, t2: t, comp: Continuous,
121 freq: NoFrequency);
122
123 const Real df
124 = std::log( x: hestonProcess_->riskFreeRate()->discount(t)
125 / hestonProcess_->riskFreeRate()->discount(t: s));
126
127 const Real eaT=std::exp(x: -a*T);
128 const Real eat=std::exp(x: -a*t);
129 const Real eas=std::exp(x: -a*s);
130 const Real iat=1.0/eat;
131 const Real ias=1.0/eas;
132
133 const Real m1 = -(dy+0.5*eta*eta)*dt - df;
134
135 const Real m2 = -rho*sigma*eta/a*(dt-1/a*eaT*(iat-ias));
136
137 const Real m3 = (r - hullWhiteProcess_->alpha(t: s))
138 *hullWhiteProcess_->B(t: s,T: t);
139
140 const Real m4 = sigma*sigma/(2*a*a)
141 *(dt + 2/a*(eat-eas) - 1/(2*a)*(eat*eat-eas*eas));
142
143 const Real m5 = -sigma*sigma/(a*a)
144 *(dt - 1/a*(1-eat*ias) - 1/(2*a)*eaT*(iat-2*ias+eat*ias*ias));
145
146 const Real mu = m1 + m2 + m3 + m4 + m5;
147
148 Array retVal(3);
149
150 const Real eta2 = hestonProcess_->sigma() * eta;
151 const Real nu
152 = hestonProcess_->kappa()*(hestonProcess_->theta() - eta*eta);
153
154 retVal[1] = x0[1] + nu*dt + eta2*std::sqrt(x: dt)
155 *(xi*dw[0]+std::sqrt(x: 1-xi*xi)*dw[1]);
156
157 if (discretization_ == BSMHullWhite) {
158 const Real v1 = eta*eta*dt
159 + sigma*sigma/(a*a)*(dt - 2/a*(1 - eat*ias)
160 + 1/(2*a)*(1 - eat*eat*ias*ias))
161 + 2*sigma*eta/a*rho*(dt - 1/a*(1 - eat*ias));
162 const Real v2 = hullWhiteProcess_->variance(t0, x0: r, dt);
163 const Real v12 = (1-eat*ias)*(sigma*eta/a*rho + sigma*sigma/(a*a))
164 - sigma*sigma/(2*a*a)*(1 - eat*eat*ias*ias);
165
166 QL_REQUIRE(v1 > 0.0 && v2 > 0.0, "zero or negative variance given");
167
168 // terminal rho must be between -maxRho and +maxRho
169 const Real rhoT
170 = std::min(a: maxRho_, b: std::max(a: -maxRho_, b: v12/std::sqrt(x: v1*v2)));
171 QL_REQUIRE( rhoT <= 1.0 && rhoT >= -1.0
172 && 1-rhoT*rhoT/(1-xi*xi) >= 0.0,
173 "invalid terminal correlation");
174
175 const Real dw_0 = dw[0];
176 const Real dw_2 = rhoT*dw[0]- rhoT*xi/std::sqrt(x: 1-xi*xi)*dw[1]
177 + std::sqrt(x: 1 - rhoT*rhoT/(1-xi*xi))*dw[2];
178
179 retVal[2] = hullWhiteProcess_->evolve(t0, x0: r, dt, dw: dw_2);
180
181 const Real vol = std::sqrt(x: v1)*dw_0;
182 retVal[0] = x0[0]*std::exp(x: mu + vol);
183 }
184 else if (discretization_ == Euler) {
185 const Real dw_2 = rho*dw[0]- rho*xi/std::sqrt(x: 1-xi*xi)*dw[1]
186 + std::sqrt(x: 1 - rho*rho/(1-xi*xi))*dw[2];
187
188 retVal[2] = hullWhiteProcess_->evolve(t0, x0: r, dt, dw: dw_2);
189
190 const Real vol = eta*std::sqrt(x: dt)*dw[0];
191 retVal[0] = x0[0]*std::exp(x: mu + vol);
192 }
193 else
194 QL_FAIL("unknown discretization scheme");
195
196 return retVal;
197 }
198
199 DiscountFactor
200 HybridHestonHullWhiteProcess::numeraire(Time t, const Array& x) const {
201
202 return hullWhiteModel_->discountBond(now: t, maturity: T_, rate: x[2]) / endDiscount_;
203 }
204
205 Real HybridHestonHullWhiteProcess::eta() const {
206 return corrEquityShortRate_;
207 }
208
209 const ext::shared_ptr<HestonProcess>&
210 HybridHestonHullWhiteProcess::hestonProcess() const {
211 return hestonProcess_;
212 }
213
214 const ext::shared_ptr<HullWhiteForwardProcess>&
215 HybridHestonHullWhiteProcess::hullWhiteProcess() const {
216 return hullWhiteProcess_;
217 }
218
219 HybridHestonHullWhiteProcess::Discretization
220 HybridHestonHullWhiteProcess::discretization() const {
221 return discretization_;
222 }
223
224 Time HybridHestonHullWhiteProcess::time(const Date& date) const {
225 return hestonProcess_->time(date);
226 }
227
228 void HybridHestonHullWhiteProcess::update() {
229 endDiscount_ = hestonProcess_->riskFreeRate()->discount(t: T_);
230 }
231}
232

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