1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2006 Banca Profilo S.p.A.
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/processes/g2process.hpp>
21#include <ql/processes/eulerdiscretization.hpp>
22
23namespace QuantLib {
24
25 G2Process::G2Process(Real a, Real sigma, Real b, Real eta, Real rho)
26 : a_(a), sigma_(sigma), b_(b), eta_(eta), rho_(rho),
27 xProcess_(new QuantLib::OrnsteinUhlenbeckProcess(a, sigma, 0.0)),
28 yProcess_(new QuantLib::OrnsteinUhlenbeckProcess(b, eta, 0.0)) {}
29
30 Size G2Process::size() const {
31 return 2;
32 }
33
34 Array G2Process::initialValues() const {
35 return { x0_, y0_ };
36 }
37
38 Array G2Process::drift(Time t, const Array& x) const {
39 return {
40 xProcess_->drift(t, x: x[0]),
41 yProcess_->drift(t, x: x[1])
42 };
43 }
44
45 Matrix G2Process::diffusion(Time, const Array&) const {
46 /* the correlation matrix is
47 | 1 rho |
48 | rho 1 |
49 whose square root (which is used here) is
50 | 1 0 |
51 | rho sqrt(1-rho^2) |
52 */
53 Matrix tmp(2,2);
54 Real sigma1 = sigma_;
55 Real sigma2 = eta_;
56 tmp[0][0] = sigma1; tmp[0][1] = 0.0;
57 tmp[1][0] = rho_*sigma1; tmp[1][1] = std::sqrt(x: 1.0-rho_*rho_)*sigma2;
58 return tmp;
59 }
60
61 Array G2Process::expectation(Time t0, const Array& x0,
62 Time dt) const {
63 return {
64 xProcess_->expectation(t0, x0: x0[0], dt),
65 yProcess_->expectation(t0, x0: x0[1], dt)
66 };
67 }
68
69 Matrix G2Process::stdDeviation(Time t0, const Array& x0, Time dt) const {
70 /* the correlation matrix is
71 | 1 rho |
72 | rho 1 |
73 whose square root (which is used here) is
74 | 1 0 |
75 | rho sqrt(1-rho^2) |
76 */
77 Matrix tmp(2,2);
78 Real sigma1 = xProcess_->stdDeviation(t: t0, x0: x0[0], dt);
79 Real sigma2 = yProcess_->stdDeviation(t: t0, x0: x0[1], dt);
80 Real expa = std::exp(x: -a_*dt), expb = std::exp(x: -b_*dt);
81 Real H = (rho_*sigma_*eta_)/(a_+b_)*(1-expa*expb);
82 Real den =
83 (0.5*sigma_*eta_)*std::sqrt(x: (1-expa*expa)*(1-expb*expb)/(a_*b_));
84 Real newRho = H/den;
85 tmp[0][0] = sigma1;
86 tmp[0][1] = 0.0;
87 tmp[1][0] = newRho*sigma2;
88 tmp[1][1] = std::sqrt(x: 1.0-newRho*newRho)*sigma2;
89 return tmp;
90 }
91
92 Matrix G2Process::covariance(Time t0, const Array& x0, Time dt) const {
93 Matrix sigma = stdDeviation(t0, x0, dt);
94 Matrix result = sigma*transpose(m: sigma);
95 return result;
96 }
97
98 Real G2Process::x0() const {
99 return x0_;
100 }
101
102 Real G2Process::y0() const {
103 return y0_;
104 }
105
106 Real G2Process::a() const {
107 return a_;
108 }
109
110 Real G2Process::sigma() const {
111 return sigma_;
112 }
113
114 Real G2Process::b() const {
115 return b_;
116 }
117
118 Real G2Process::eta() const {
119 return eta_;
120 }
121
122 Real G2Process::rho() const {
123 return rho_;
124 }
125
126
127 G2ForwardProcess::G2ForwardProcess(Real a, Real sigma, Real b, Real eta, Real rho)
128 : a_(a), sigma_(sigma), b_(b), eta_(eta), rho_(rho),
129 xProcess_(new QuantLib::OrnsteinUhlenbeckProcess(a, sigma, 0.0)),
130 yProcess_(new QuantLib::OrnsteinUhlenbeckProcess(b, eta, 0.0)) {}
131
132 Size G2ForwardProcess::size() const {
133 return 2;
134 }
135
136 Array G2ForwardProcess::initialValues() const {
137 return { x0_, y0_ };
138 }
139
140 Array G2ForwardProcess::drift(Time t, const Array& x) const {
141 return {
142 xProcess_->drift(t, x: x[0]) + xForwardDrift(t, T: T_),
143 yProcess_->drift(t, x: x[1]) + yForwardDrift(t, T: T_)
144 };
145 }
146
147 Matrix G2ForwardProcess::diffusion(Time, const Array&) const {
148 Matrix tmp(2,2);
149 Real sigma1 = sigma_;
150 Real sigma2 = eta_;
151 tmp[0][0] = sigma1; tmp[0][1] = 0.0;
152 tmp[1][0] = rho_*sigma1; tmp[1][1] = std::sqrt(x: 1.0-rho_*rho_)*sigma2;
153 return tmp;
154 }
155
156 Array G2ForwardProcess::expectation(Time t0, const Array& x0,
157 Time dt) const {
158 return {
159 xProcess_->expectation(t0, x0: x0[0], dt) - Mx_T(s: t0, t: t0+dt, T: T_),
160 yProcess_->expectation(t0, x0: x0[1], dt) - My_T(s: t0, t: t0+dt, T: T_)
161 };
162 }
163
164 Matrix G2ForwardProcess::stdDeviation(Time t0, const Array& x0, Time dt) const {
165 Matrix tmp(2,2);
166 Real sigma1 = xProcess_->stdDeviation(t: t0, x0: x0[0], dt);
167 Real sigma2 = yProcess_->stdDeviation(t: t0, x0: x0[1], dt);
168 Real expa = std::exp(x: -a_*dt), expb = std::exp(x: -b_*dt);
169 Real H = (rho_*sigma_*eta_)/(a_+b_)*(1-expa*expb);
170 Real den =
171 (0.5*sigma_*eta_)*std::sqrt(x: (1-expa*expa)*(1-expb*expb)/(a_*b_));
172 Real newRho = H/den;
173 tmp[0][0] = sigma1;
174 tmp[0][1] = 0.0;
175 tmp[1][0] = newRho*sigma2;
176 tmp[1][1] = std::sqrt(x: 1.0-newRho*newRho)*sigma2;
177 return tmp;
178 }
179
180 Matrix G2ForwardProcess::covariance(Time t0, const Array& x0, Time dt) const {
181 Matrix sigma = stdDeviation(t0, x0, dt);
182 Matrix result = sigma*transpose(m: sigma);
183 return result;
184 }
185
186 Real G2ForwardProcess::xForwardDrift(Time t, Time T) const {
187 Real expatT = std::exp(x: -a_*(T-t));
188 Real expbtT = std::exp(x: -b_*(T-t));
189
190 return -(sigma_*sigma_/a_) * (1-expatT)
191 - (rho_*sigma_*eta_/b_) * (1-expbtT);
192 }
193
194 Real G2ForwardProcess::yForwardDrift(Time t, Time T) const {
195 Real expatT = std::exp(x: -a_*(T-t));
196 Real expbtT = std::exp(x: -b_*(T-t));
197
198 return -(eta_*eta_/b_) * (1-expbtT)
199 - (rho_*sigma_*eta_/a_) * (1-expatT);
200 }
201
202 Real G2ForwardProcess::Mx_T(Real s, Real t, Real T) const {
203 Real M;
204 M = ( (sigma_*sigma_)/(a_*a_) + (rho_*sigma_*eta_)/(a_*b_) )
205 * (1-std::exp(x: -a_*(t-s)));
206 M += -(sigma_*sigma_)/(2*a_*a_) *
207 (std::exp(x: -a_*(T-t))-std::exp(x: -a_*(T+t-2*s)));
208 M += -(rho_*sigma_*eta_)/(b_*(a_+b_))
209 * (std::exp(x: -b_*(T-t)) -std::exp(x: -b_*T-a_*t+(a_+b_)*s));
210 return M;
211 }
212
213 Real G2ForwardProcess::My_T(Real s, Real t, Real T) const {
214 Real M;
215 M = ( (eta_*eta_)/(b_*b_) + (rho_*sigma_*eta_)/(a_*b_) )
216 * (1-std::exp(x: -b_*(t-s)));
217 M += -(eta_*eta_)/(2*b_*b_) *
218 (std::exp(x: -b_*(T-t))-std::exp(x: -b_*(T+t-2*s)));
219 M += -(rho_*sigma_*eta_)/(a_*(a_+b_))
220 * (std::exp(x: -a_*(T-t))-std::exp(x: -a_*T-b_*t+(a_+b_)*s));
221 return M;
222 }
223
224}
225
226

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