11328508063

Committed 14 Oct 2024 01:24PM UTC coverage: 72.674% (+0.05%) from 72.621%

Build # 11328508063

Build Type

push

github

Committed by

lballabio

Commit Message

Pin Ubuntu version in coverage workflow

Run Details

55158 of 75898 relevant lines covered (72.67%)

8693210.46 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

0.0

/ql/processes/jointstochasticprocess.cpp

/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */

/*
 Copyright (C) 2007 Klaus Spanderen

 This file is part of QuantLib, a free-software/open-source library
 for financial quantitative analysts and developers - http://quantlib.org/

 QuantLib is free software: you can redistribute it and/or modify it
 under the terms of the QuantLib license.  You should have received a
 copy of the license along with this program; if not, please email
 <quantlib-dev@lists.sf.net>. The license is also available online at
 <http://quantlib.org/license.shtml>.

 This program is distributed in the hope that it will be useful, but WITHOUT
 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 FOR A PARTICULAR PURPOSE.  See the license for more details.
*/

/*! \file jointstochasticprocess.cpp
    \brief multi model process for hybrid products
*/

#include <ql/math/matrixutilities/pseudosqrt.hpp>
#include <ql/math/matrixutilities/svd.hpp>
#include <ql/processes/jointstochasticprocess.hpp>
#include <utility>

namespace QuantLib {

    JointStochasticProcess::JointStochasticProcess(
        std::vector<ext::shared_ptr<StochasticProcess> > l, Size factors)
    : l_(std::move(l)), factors_(factors) {

        for (const auto& iter : l_) {
            registerWith(iter);
        }

        vsize_.reserve   (l_.size()+1);
        vfactors_.reserve(l_.size()+1);

        for (const auto& iter : l_) {
            vsize_.push_back(size_);
            size_ += iter->size();

            vfactors_.push_back(modelFactors_);
            modelFactors_ += iter->factors();
        }

        vsize_.push_back(size_);
        vfactors_.push_back(modelFactors_);

        if (factors_ == Null<Size>()) {
            factors_ = modelFactors_;
        } else {
            QL_REQUIRE(factors_ <= size_, "too many factors given");
        }
    }

    Size JointStochasticProcess::size() const {
        return size_;
    }

    Size JointStochasticProcess::factors() const {
        return factors_;
    }

    Array JointStochasticProcess::slice(const Array& x,
                                        Size i) const {
        // cut out the ith process' variables
        Size n = vsize_[i+1]-vsize_[i];
        Array y(n);
        std::copy(x.begin()+vsize_[i], x.begin()+vsize_[i+1], y.begin());
        return y;
    }

    Array JointStochasticProcess::initialValues() const {
        Array retVal(size());

        for (auto iter = l_.begin(); iter != l_.end(); ++iter) {
            const Array& pInitValues = (*iter)->initialValues();

            std::copy(pInitValues.begin(), pInitValues.end(),
                      retVal.begin()+vsize_[iter - l_.begin()]);
        }

        return retVal;
    }


    Array JointStochasticProcess::drift(Time t,
                                        const Array& x) const {
        Array retVal(size());

        for (Size i=0; i < l_.size(); ++i) {

            const Array& pDrift = l_[i]->drift(t, slice(x,i));

            std::copy(pDrift.begin(), pDrift.end(),
                      retVal.begin()+vsize_[i]);
        }

        return retVal;
    }

    Array JointStochasticProcess::expectation(Time t0,
                                              const Array& x0,
                                              Time dt) const {
        Array retVal(size());

        for (Size i=0; i < l_.size(); ++i) {

            const Array& pExpectation = l_[i]->expectation(t0, slice(x0,i), dt);

            std::copy(pExpectation.begin(), pExpectation.end(),
                      retVal.begin()+ vsize_[i]);
        }

        return retVal;
    }


    Matrix JointStochasticProcess::diffusion(Time t, const Array& x) const {
        // might need some improvement in the future
        const Time dt = 0.001;
        return pseudoSqrt(covariance(t, x, dt)/dt);
    }


    Matrix JointStochasticProcess::covariance(Time t0,
                                              const Array& x0,
                                              Time dt) const {

        // get the model intrinsic covariance matrix
        Matrix retVal(size(), size(), 0.0);

        for (Size j=0; j < l_.size(); ++j) {

            const Size vs = vsize_[j];
            const Matrix& pCov = l_[j]->covariance(t0, slice(x0,j), dt);

            for (Size i=0; i < pCov.rows(); ++i) {
                std::copy(pCov.row_begin(i), pCov.row_end(i),
                          retVal.row_begin(vs+i) + vs);
            }
        }

        // add the cross model covariance matrix
        const Array& volatility = Sqrt(retVal.diagonal());
        Matrix crossModelCovar = this->crossModelCorrelation(t0, x0);

        for (Size i=0; i < size(); ++i) {
            for (Size j=0; j < size(); ++j) {
                crossModelCovar[i][j] *= volatility[i]*volatility[j];
            }
        }

        retVal += crossModelCovar;

        return retVal;
    }


    Matrix JointStochasticProcess::stdDeviation(Time t0,
                                                const Array& x0,
                                                Time dt) const {
        return pseudoSqrt(covariance(t0, x0, dt));
    }


    Array JointStochasticProcess::apply(const Array& x0,
                                        const Array& dx) const {
        Array retVal(size());

        for (Size i=0; i < l_.size(); ++i) {
            const Array& pApply = l_[i]->apply(slice(x0,i), slice(dx,i));

            std::copy(pApply.begin(), pApply.end(),
                      retVal.begin()+vsize_[i]);
        }

        return retVal;
    }

    Array JointStochasticProcess::evolve(
        Time t0, const Array& x0, Time dt, const Array& dw) const {
        Array dv(modelFactors_);

        if (   correlationIsStateDependent()
            || correlationCache_.count(CachingKey(t0, dt)) == 0) {
            Matrix cov  = covariance(t0, x0, dt);

            const Array& sqrtDiag = Sqrt(cov.diagonal());
            for (Size i=0; i < cov.rows(); ++i) {
                for (Size j=i; j < cov.columns(); ++j) {
                    const Real div = sqrtDiag[i]*sqrtDiag[j];

                    cov[i][j] = cov[j][i] = ( div > 0) ? Real(cov[i][j]/div) : 0.0;
                }
            }

            Matrix diff(size(), modelFactors_, 0.0);

            for (Size j = 0; j < l_.size(); ++j) {
                const Size vs = vsize_   [j];
                const Size vf = vfactors_[j];

                Matrix stdDev = l_[j]->stdDeviation(t0, slice(x0,j), dt);

                for (Size i=0; i < stdDev.rows(); ++i) {
                    const Volatility vol = std::sqrt(
                        std::inner_product(stdDev.row_begin(i),
                                           stdDev.row_end(i),
                                           stdDev.row_begin(i), Real(0.0)));
                    if (vol > 0.0) {
                        std::transform(stdDev.row_begin(i), stdDev.row_end(i),
                                       stdDev.row_begin(i),
                                       [=](Real x) -> Real { return x / vol; });
                    }
                    else {
                        // keep the svd happy
                        std::fill(stdDev.row_begin(i), stdDev.row_end(i),
                                  100*i*QL_EPSILON);
                    }
                }

                SVD svd(stdDev);
                const Array& s = svd.singularValues();
                Matrix w(s.size(), s.size(), 0.0);
                for (Size i=0; i < s.size(); ++i) {
                    if (std::fabs(s[i]) > std::sqrt(QL_EPSILON)) {
                        w[i][i] = 1.0/s[i];
                    }
                }

                const Matrix inv = svd.U() * w * transpose(svd.V());

                for (Size i=0; i < stdDev.rows(); ++i) {
                    std::copy(inv.row_begin(i), inv.row_end(i),
                              diff.row_begin(i+vs)+vf);
                }
            }

            Matrix rs = rankReducedSqrt(cov, factors_, 1.0,
                                        SalvagingAlgorithm::Spectral);

            if (rs.columns() < factors_) {
                // less eigenvalues than expected factors.
                // fill the rest with zero's.
                Matrix tmp = Matrix(cov.rows(), factors_, 0.0);
                for (Size i=0; i < cov.rows(); ++i) {
                    std::copy(rs.row_begin(i), rs.row_end(i),
                              tmp.row_begin(i));
                }
                rs = tmp;
            }

            const Matrix m = transpose(diff) * rs;

            if (!correlationIsStateDependent()) {
                correlationCache_[CachingKey(t0,dt)] = m;
            }
            dv = m*dw;
        }
        else {
            if (!correlationIsStateDependent()) {
                dv = correlationCache_[CachingKey(t0,dt)] * dw;
            }
        }

        this->preEvolve(t0, x0, dt, dv);


        Array retVal(size());
        for (auto iter = l_.begin(); iter != l_.end(); ++iter) {
            const Size i = iter - l_.begin();

            Array dz((*iter)->factors());
            std::copy(dv.begin()+vfactors_[i],
                      dv.begin()+vfactors_[i] + (*iter)->factors(),
                      dz.begin());
            Array x((*iter)->size());
            std::copy(x0.begin()+vsize_[i],
                      x0.begin()+vsize_[i] + (*iter)->size(),
                      x.begin());
            const Array r = (*iter)->evolve(t0, x, dt, dz);
            std::copy(r.begin(), r.end(), retVal.begin()+vsize_[i]);
        }

        return this->postEvolve(t0, x0, dt, dv, retVal);
    }

    const std::vector<ext::shared_ptr<StochasticProcess> > &
                          JointStochasticProcess::constituents() const {
        return l_;
    }

    Time JointStochasticProcess::time(const Date& date) const {
        QL_REQUIRE(!l_.empty(), "process list is empty");

        return l_[0]->time(date);
    }

    void JointStochasticProcess::update() {
        // clear all caches
        correlationCache_.clear();

        this->StochasticProcess::update();
    }
}

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
29	namespace 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(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>()) {	×
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(x.begin()+vsize_[i], x.begin()+vsize_[i+1], 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(pInitValues.begin(), pInitValues.end(),	×
84	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, slice(x,i));	×
98
99	std::copy(pDrift.begin(), pDrift.end(),	×
100	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, 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
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(t, 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, 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
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(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
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(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	100iQL_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,	×
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(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
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, 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> > &
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	}

lballabio / QuantLib / 11328508063

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous