QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
jointstochasticprocess.cpp
Go to the documentation of this file.
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
27#include <utility>
28
29namespace QuantLib {
30
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(iter);
37 }
38
39 vsize_.reserve (l_.size()+1);
40 vfactors_.reserve(l_.size()+1);
41
42 for (const auto& iter : l_) {
43 vsize_.push_back(size_);
44 size_ += iter->size();
45
46 vfactors_.push_back(modelFactors_);
47 modelFactors_ += iter->factors();
48 }
49
50 vsize_.push_back(size_);
51 vfactors_.push_back(modelFactors_);
52
53 if (factors_ == Null<Size>()) {
55 } else {
56 QL_REQUIRE(factors_ <= size_, "too many factors given");
57 }
58 }
59
61 return size_;
62 }
63
65 return factors_;
66 }
67
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(x.begin()+vsize_[i], x.begin()+vsize_[i+1], y.begin());
74 return y;
75 }
76
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(pInitValues.begin(), pInitValues.end(),
84 retVal.begin()+vsize_[iter - l_.begin()]);
85 }
86
87 return retVal;
88 }
89
90
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, slice(x,i));
98
99 std::copy(pDrift.begin(), pDrift.end(),
100 retVal.begin()+vsize_[i]);
101 }
102
103 return retVal;
104 }
105
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, slice(x0,i), dt);
114
115 std::copy(pExpectation.begin(), pExpectation.end(),
116 retVal.begin()+ vsize_[i]);
117 }
118
119 return retVal;
120 }
121
122
124 // might need some improvement in the future
125 const Time dt = 0.001;
126 return pseudoSqrt(covariance(t, x, dt)/dt);
127 }
128
129
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, slice(x0,j), dt);
141
142 for (Size i=0; i < pCov.rows(); ++i) {
143 std::copy(pCov.row_begin(i), pCov.row_end(i),
144 retVal.row_begin(vs+i) + vs);
145 }
146 }
147
148 // add the cross model covariance matrix
149 const Array& volatility = Sqrt(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
165 const Array& x0,
166 Time dt) const {
167 return pseudoSqrt(covariance(t0, x0, dt));
168 }
169
170
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(slice(x0,i), slice(dx,i));
177
178 std::copy(pApply.begin(), pApply.end(),
179 retVal.begin()+vsize_[i]);
180 }
181
182 return retVal;
183 }
184
186 Time t0, const Array& x0, Time dt, const Array& dw) const {
188
190 || correlationCache_.count(CachingKey(t0, dt)) == 0) {
191 Matrix cov = covariance(t0, x0, dt);
192
193 const Array& sqrtDiag = Sqrt(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, slice(x0,j), dt);
209
210 for (Size i=0; i < stdDev.rows(); ++i) {
211 const Volatility vol = std::sqrt(
212 std::inner_product(stdDev.row_begin(i),
213 stdDev.row_end(i),
214 stdDev.row_begin(i), Real(0.0)));
215 if (vol > 0.0) {
216 std::transform(stdDev.row_begin(i), stdDev.row_end(i),
217 stdDev.row_begin(i),
218 [=](Real x) -> Real { return x / vol; });
219 }
220 else {
221 // keep the svd happy
222 std::fill(stdDev.row_begin(i), stdDev.row_end(i),
223 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(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(svd.V());
237
238 for (Size i=0; i < stdDev.rows(); ++i) {
239 std::copy(inv.row_begin(i), inv.row_end(i),
240 diff.row_begin(i+vs)+vf);
241 }
242 }
243
244 Matrix rs = rankReducedSqrt(cov, factors_, 1.0,
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(rs.row_begin(i), rs.row_end(i),
253 tmp.row_begin(i));
254 }
255 rs = tmp;
256 }
257
258 const Matrix m = transpose(diff) * rs;
259
261 correlationCache_[CachingKey(t0,dt)] = m;
262 }
263 dv = m*dw;
264 }
265 else {
267 dv = correlationCache_[CachingKey(t0,dt)] * dw;
268 }
269 }
270
271 this->preEvolve(t0, x0, dt, 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(dv.begin()+vfactors_[i],
280 dv.begin()+vfactors_[i] + (*iter)->factors(),
281 dz.begin());
282 Array x((*iter)->size());
283 std::copy(x0.begin()+vsize_[i],
284 x0.begin()+vsize_[i] + (*iter)->size(),
285 x.begin());
286 const Array r = (*iter)->evolve(t0, x, dt, dz);
287 std::copy(r.begin(), r.end(), retVal.begin()+vsize_[i]);
288 }
289
290 return this->postEvolve(t0, x0, dt, dv, retVal);
291 }
292
293 const std::vector<ext::shared_ptr<StochasticProcess> > &
295 return l_;
296 }
297
299 QL_REQUIRE(!l_.empty(), "process list is empty");
300
301 return l_[0]->time(date);
302 }
303
305 // clear all caches
306 correlationCache_.clear();
307
309 }
310}
1-D array used in linear algebra.
Definition: array.hpp:52
const_iterator end() const
Definition: array.hpp:511
const_iterator begin() const
Definition: array.hpp:503
Concrete date class.
Definition: date.hpp:125
std::vector< ext::shared_ptr< StochasticProcess > > l_
Array slice(const Array &x, Size i) const
virtual Matrix crossModelCorrelation(Time t0, const Array &x0) const =0
Array drift(Time t, const Array &x) const override
returns the drift part of the equation, i.e.,
Size size() const override
returns the number of dimensions of the stochastic process
Array evolve(Time t0, const Array &x0, Time dt, const Array &dw) const override
Matrix diffusion(Time t, const Array &x) const override
returns the diffusion part of the equation, i.e.
Time time(const Date &date) const override
Matrix stdDeviation(Time t0, const Array &x0, Time dt) const override
std::map< CachingKey, Matrix > correlationCache_
virtual Array postEvolve(Time t0, const Array &x0, Time dt, const Array &dw, const Array &y0) const =0
virtual void preEvolve(Time t0, const Array &x0, Time dt, const Array &dw) const =0
virtual bool correlationIsStateDependent() const =0
Size factors() const override
returns the number of independent factors of the process
Matrix covariance(Time t0, const Array &x0, Time dt) const override
Array expectation(Time t0, const Array &x0, Time dt) const override
Array initialValues() const override
returns the initial values of the state variables
Array apply(const Array &x0, const Array &dx) const override
JointStochasticProcess(std::vector< ext::shared_ptr< StochasticProcess > > l, Size factors=Null< Size >())
const std::vector< ext::shared_ptr< StochasticProcess > > & constituents() const
Matrix used in linear algebra.
Definition: matrix.hpp:41
const_row_iterator row_begin(Size i) const
Definition: matrix.hpp:360
Size rows() const
Definition: matrix.hpp:504
Array diagonal() const
Definition: matrix.hpp:492
Size columns() const
Definition: matrix.hpp:508
const_row_iterator row_end(Size i) const
Definition: matrix.hpp:378
template class providing a null value for a given type.
Definition: null.hpp:76
std::pair< iterator, bool > registerWith(const ext::shared_ptr< Observable > &)
Definition: observable.hpp:228
Singular value decomposition.
Definition: svd.hpp:54
const Matrix & V() const
Definition: svd.cpp:489
const Array & singularValues() const
Definition: svd.cpp:493
const Matrix & U() const
Definition: svd.cpp:485
const DefaultType & t
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
Definition: errors.hpp:117
#define QL_EPSILON
Definition: qldefines.hpp:178
Real Time
continuous quantity with 1-year units
Definition: types.hpp:62
QL_REAL Real
real number
Definition: types.hpp:50
Real Volatility
volatility
Definition: types.hpp:78
std::size_t Size
size of a container
Definition: types.hpp:58
multi model process for hybrid products
Definition: any.hpp:35
Matrix pseudoSqrt(const Matrix &matrix, SalvagingAlgorithm::Type sa)
Definition: pseudosqrt.cpp:347
Matrix rankReducedSqrt(const Matrix &matrix, Size maxRank, Real componentRetainedPercentage, SalvagingAlgorithm::Type sa)
Definition: pseudosqrt.cpp:427
Array Sqrt(const Array &v)
Definition: array.hpp:847
Matrix transpose(const Matrix &m)
Definition: matrix.hpp:700
STL namespace.
ext::shared_ptr< YieldTermStructure > r
pseudo square root of a real symmetric matrix
singular value decomposition