1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2007 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 jointstochasticprocess.cpp
21 \brief multi model process for hybrid products
22*/
23
24#include <ql/math/matrixutilities/pseudosqrt.hpp>
25#include <ql/math/matrixutilities/svd.hpp>
26#include <ql/processes/jointstochasticprocess.hpp>
27#include <utility>
28
29namespace QuantLib {
30
31 JointStochasticProcess::JointStochasticProcess(
32 std::vector<ext::shared_ptr<StochasticProcess> > l, Size factors)
33 : l_(std::move(l)), factors_(factors) {
34
35 for (const auto& iter : l_) {
36 registerWith(h: iter);
37 }
38
39 vsize_.reserve (n: l_.size()+1);
40 vfactors_.reserve(n: l_.size()+1);
41
42 for (const auto& iter : l_) {
43 vsize_.push_back(x: size_);
44 size_ += iter->size();
45
46 vfactors_.push_back(x: modelFactors_);
47 modelFactors_ += iter->factors();
48 }
49
50 vsize_.push_back(x: size_);
51 vfactors_.push_back(x: modelFactors_);
52
53 if (factors_ == Null<Size>()) {
54 factors_ = modelFactors_;
55 } else {
56 QL_REQUIRE(factors_ <= size_, "too many factors given");
57 }
58 }
59
60 Size JointStochasticProcess::size() const {
61 return size_;
62 }
63
64 Size JointStochasticProcess::factors() const {
65 return factors_;
66 }
67
68 Array JointStochasticProcess::slice(const Array& x,
69 Size i) const {
70 // cut out the ith process' variables
71 Size n = vsize_[i+1]-vsize_[i];
72 Array y(n);
73 std::copy(first: x.begin()+vsize_[i], last: x.begin()+vsize_[i+1], result: y.begin());
74 return y;
75 }
76
77 Array JointStochasticProcess::initialValues() const {
78 Array retVal(size());
79
80 for (auto iter = l_.begin(); iter != l_.end(); ++iter) {
81 const Array& pInitValues = (*iter)->initialValues();
82
83 std::copy(first: pInitValues.begin(), last: pInitValues.end(),
84 result: retVal.begin()+vsize_[iter - l_.begin()]);
85 }
86
87 return retVal;
88 }
89
90
91 Array JointStochasticProcess::drift(Time t,
92 const Array& x) const {
93 Array retVal(size());
94
95 for (Size i=0; i < l_.size(); ++i) {
96
97 const Array& pDrift = l_[i]->drift(t, x: slice(x,i));
98
99 std::copy(first: pDrift.begin(), last: pDrift.end(),
100 result: retVal.begin()+vsize_[i]);
101 }
102
103 return retVal;
104 }
105
106 Array JointStochasticProcess::expectation(Time t0,
107 const Array& x0,
108 Time dt) const {
109 Array retVal(size());
110
111 for (Size i=0; i < l_.size(); ++i) {
112
113 const Array& pExpectation = l_[i]->expectation(t0, x0: slice(x: x0,i), dt);
114
115 std::copy(first: pExpectation.begin(), last: pExpectation.end(),
116 result: retVal.begin()+ vsize_[i]);
117 }
118
119 return retVal;
120 }
121
122
123 Matrix JointStochasticProcess::diffusion(Time t, const Array& x) const {
124 // might need some improvement in the future
125 const Time dt = 0.001;
126 return pseudoSqrt(covariance(t0: t, x0: x, dt)/dt);
127 }
128
129
130 Matrix JointStochasticProcess::covariance(Time t0,
131 const Array& x0,
132 Time dt) const {
133
134 // get the model intrinsic covariance matrix
135 Matrix retVal(size(), size(), 0.0);
136
137 for (Size j=0; j < l_.size(); ++j) {
138
139 const Size vs = vsize_[j];
140 const Matrix& pCov = l_[j]->covariance(t0, x0: slice(x: x0,i: j), dt);
141
142 for (Size i=0; i < pCov.rows(); ++i) {
143 std::copy(first: pCov.row_begin(i), last: pCov.row_end(i),
144 result: retVal.row_begin(i: vs+i) + vs);
145 }
146 }
147
148 // add the cross model covariance matrix
149 const Array& volatility = Sqrt(v: retVal.diagonal());
150 Matrix crossModelCovar = this->crossModelCorrelation(t0, x0);
151
152 for (Size i=0; i < size(); ++i) {
153 for (Size j=0; j < size(); ++j) {
154 crossModelCovar[i][j] *= volatility[i]*volatility[j];
155 }
156 }
157
158 retVal += crossModelCovar;
159
160 return retVal;
161 }
162
163
164 Matrix JointStochasticProcess::stdDeviation(Time t0,
165 const Array& x0,
166 Time dt) const {
167 return pseudoSqrt(covariance(t0, x0, dt));
168 }
169
170
171 Array JointStochasticProcess::apply(const Array& x0,
172 const Array& dx) const {
173 Array retVal(size());
174
175 for (Size i=0; i < l_.size(); ++i) {
176 const Array& pApply = l_[i]->apply(x0: slice(x: x0,i), dx: slice(x: dx,i));
177
178 std::copy(first: pApply.begin(), last: pApply.end(),
179 result: retVal.begin()+vsize_[i]);
180 }
181
182 return retVal;
183 }
184
185 Array JointStochasticProcess::evolve(
186 Time t0, const Array& x0, Time dt, const Array& dw) const {
187 Array dv(modelFactors_);
188
189 if ( correlationIsStateDependent()
190 || correlationCache_.count(x: CachingKey(t0, dt)) == 0) {
191 Matrix cov = covariance(t0, x0, dt);
192
193 const Array& sqrtDiag = Sqrt(v: cov.diagonal());
194 for (Size i=0; i < cov.rows(); ++i) {
195 for (Size j=i; j < cov.columns(); ++j) {
196 const Real div = sqrtDiag[i]*sqrtDiag[j];
197
198 cov[i][j] = cov[j][i] = ( div > 0) ? Real(cov[i][j]/div) : 0.0;
199 }
200 }
201
202 Matrix diff(size(), modelFactors_, 0.0);
203
204 for (Size j = 0; j < l_.size(); ++j) {
205 const Size vs = vsize_ [j];
206 const Size vf = vfactors_[j];
207
208 Matrix stdDev = l_[j]->stdDeviation(t0, x0: slice(x: x0,i: j), dt);
209
210 for (Size i=0; i < stdDev.rows(); ++i) {
211 const Volatility vol = std::sqrt(
212 x: std::inner_product(first1: stdDev.row_begin(i),
213 last1: stdDev.row_end(i),
214 first2: stdDev.row_begin(i), init: Real(0.0)));
215 if (vol > 0.0) {
216 std::transform(first: stdDev.row_begin(i), last: stdDev.row_end(i),
217 result: stdDev.row_begin(i),
218 unary_op: [=](Real x) -> Real { return x / vol; });
219 }
220 else {
221 // keep the svd happy
222 std::fill(first: stdDev.row_begin(i), last: stdDev.row_end(i),
223 value: 100*i*QL_EPSILON);
224 }
225 }
226
227 SVD svd(stdDev);
228 const Array& s = svd.singularValues();
229 Matrix w(s.size(), s.size(), 0.0);
230 for (Size i=0; i < s.size(); ++i) {
231 if (std::fabs(x: s[i]) > std::sqrt(QL_EPSILON)) {
232 w[i][i] = 1.0/s[i];
233 }
234 }
235
236 const Matrix inv = svd.U() * w * transpose(m: svd.V());
237
238 for (Size i=0; i < stdDev.rows(); ++i) {
239 std::copy(first: inv.row_begin(i), last: inv.row_end(i),
240 result: diff.row_begin(i: i+vs)+vf);
241 }
242 }
243
244 Matrix rs = rankReducedSqrt(cov, maxRank: factors_, componentRetainedPercentage: 1.0,
245 SalvagingAlgorithm::Spectral);
246
247 if (rs.columns() < factors_) {
248 // less eigenvalues than expected factors.
249 // fill the rest with zero's.
250 Matrix tmp = Matrix(cov.rows(), factors_, 0.0);
251 for (Size i=0; i < cov.rows(); ++i) {
252 std::copy(first: rs.row_begin(i), last: rs.row_end(i),
253 result: tmp.row_begin(i));
254 }
255 rs = tmp;
256 }
257
258 const Matrix m = transpose(m: diff) * rs;
259
260 if (!correlationIsStateDependent()) {
261 correlationCache_[CachingKey(t0,dt)] = m;
262 }
263 dv = m*dw;
264 }
265 else {
266 if (!correlationIsStateDependent()) {
267 dv = correlationCache_[CachingKey(t0,dt)] * dw;
268 }
269 }
270
271 this->preEvolve(t0, x0, dt, dw: dv);
272
273
274 Array retVal(size());
275 for (auto iter = l_.begin(); iter != l_.end(); ++iter) {
276 const Size i = iter - l_.begin();
277
278 Array dz((*iter)->factors());
279 std::copy(first: dv.begin()+vfactors_[i],
280 last: dv.begin()+vfactors_[i] + (*iter)->factors(),
281 result: dz.begin());
282 Array x((*iter)->size());
283 std::copy(first: x0.begin()+vsize_[i],
284 last: x0.begin()+vsize_[i] + (*iter)->size(),
285 result: x.begin());
286 const Array r = (*iter)->evolve(t0, x0: x, dt, dw: dz);
287 std::copy(first: r.begin(), last: r.end(), result: retVal.begin()+vsize_[i]);
288 }
289
290 return this->postEvolve(t0, x0, dt, dw: dv, y0: retVal);
291 }
292
293 const std::vector<ext::shared_ptr<StochasticProcess> > &
294 JointStochasticProcess::constituents() const {
295 return l_;
296 }
297
298 Time JointStochasticProcess::time(const Date& date) const {
299 QL_REQUIRE(!l_.empty(), "process list is empty");
300
301 return l_[0]->time(date);
302 }
303
304 void JointStochasticProcess::update() {
305 // clear all caches
306 correlationCache_.clear();
307
308 this->StochasticProcess::update();
309 }
310}
311

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