/src/quantlib/ql/processes/g2process.cpp

Source (jump to first uncovered line)
/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */

/*
 Copyright (C) 2006 Banca Profilo S.p.A.

 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.
*/

#include <ql/processes/g2process.hpp>
#include <ql/processes/eulerdiscretization.hpp>

namespace QuantLib {

    G2Process::G2Process(Real a, Real sigma, Real b, Real eta, Real rho)
    : a_(a), sigma_(sigma), b_(b), eta_(eta), rho_(rho),
      xProcess_(new QuantLib::OrnsteinUhlenbeckProcess(a, sigma, 0.0)),
      yProcess_(new QuantLib::OrnsteinUhlenbeckProcess(b, eta, 0.0)) {}

    Size G2Process::size() const {
        return 2;
    }

    Array G2Process::initialValues() const {
        return { x0_, y0_ };
    }

    Array G2Process::drift(Time t, const Array& x) const {
        return {
            xProcess_->drift(t, x[0]),
            yProcess_->drift(t, x[1])
        };
    }

    Matrix G2Process::diffusion(Time, const Array&) const {
        /* the correlation matrix is
           |  1   rho |
           | rho   1  |
           whose square root (which is used here) is
           |  1          0       |
           | rho   sqrt(1-rho^2) |
        */
        Matrix tmp(2,2);
        Real sigma1 = sigma_;
        Real sigma2 = eta_;
        tmp[0][0] = sigma1;       tmp[0][1] = 0.0;
        tmp[1][0] = rho_*sigma1;  tmp[1][1] = std::sqrt(1.0-rho_*rho_)*sigma2;
        return tmp;
    }

    Array G2Process::expectation(Time t0, const Array& x0,
                                 Time dt) const {
        return {
            xProcess_->expectation(t0, x0[0], dt),
            yProcess_->expectation(t0, x0[1], dt)
        };
    }

    Matrix G2Process::stdDeviation(Time t0, const Array& x0, Time dt) const {
        /* the correlation matrix is
           |  1   rho |
           | rho   1  |
           whose square root (which is used here) is
           |  1          0       |
           | rho   sqrt(1-rho^2) |
        */
        Matrix tmp(2,2);
        Real sigma1 = xProcess_->stdDeviation(t0, x0[0], dt);
        Real sigma2 = yProcess_->stdDeviation(t0, x0[1], dt);
        Real expa = std::exp(-a_*dt), expb = std::exp(-b_*dt);
        Real H = (rho_*sigma_*eta_)/(a_+b_)*(1-expa*expb);
        Real den =
            (0.5*sigma_*eta_)*std::sqrt((1-expa*expa)*(1-expb*expb)/(a_*b_));
        Real newRho = H/den;
        tmp[0][0] = sigma1;
        tmp[0][1] = 0.0;
        tmp[1][0] = newRho*sigma2;
        tmp[1][1] = std::sqrt(1.0-newRho*newRho)*sigma2;
        return tmp;
    }

    Matrix G2Process::covariance(Time t0, const Array& x0, Time dt) const {
        Matrix sigma = stdDeviation(t0, x0, dt);
        Matrix result = sigma*transpose(sigma);
        return result;
    }

    Real G2Process::x0() const {
        return x0_;
    }

    Real G2Process::y0() const {
        return y0_;
    }

    Real G2Process::a() const {
        return a_;
    }

    Real G2Process::sigma() const {
        return sigma_;
    }

    Real G2Process::b() const {
        return b_;
    }

    Real G2Process::eta() const {
        return eta_;
    }

    Real G2Process::rho() const {
        return rho_;
    }


    G2ForwardProcess::G2ForwardProcess(Real a, Real sigma, Real b, Real eta, Real rho)
    : a_(a), sigma_(sigma), b_(b), eta_(eta), rho_(rho),
      xProcess_(new QuantLib::OrnsteinUhlenbeckProcess(a, sigma, 0.0)),
      yProcess_(new QuantLib::OrnsteinUhlenbeckProcess(b, eta, 0.0)) {}

    Size G2ForwardProcess::size() const {
        return 2;
    }

    Array G2ForwardProcess::initialValues() const {
        return { x0_, y0_ };
    }

    Array G2ForwardProcess::drift(Time t, const Array& x) const {
        return {
            xProcess_->drift(t, x[0]) + xForwardDrift(t, T_),
            yProcess_->drift(t, x[1]) + yForwardDrift(t, T_)
        };
    }

    Matrix G2ForwardProcess::diffusion(Time, const Array&) const {
        Matrix tmp(2,2);
        Real sigma1 = sigma_;
        Real sigma2 = eta_;
        tmp[0][0] = sigma1;       tmp[0][1] = 0.0;
        tmp[1][0] = rho_*sigma1;  tmp[1][1] = std::sqrt(1.0-rho_*rho_)*sigma2;
        return tmp;
    }

    Array G2ForwardProcess::expectation(Time t0, const Array& x0,
                                        Time dt) const {
        return {
            xProcess_->expectation(t0, x0[0], dt) - Mx_T(t0, t0+dt, T_),
            yProcess_->expectation(t0, x0[1], dt) - My_T(t0, t0+dt, T_)
        };
    }

    Matrix G2ForwardProcess::stdDeviation(Time t0, const Array& x0, Time dt) const {
        Matrix tmp(2,2);
        Real sigma1 = xProcess_->stdDeviation(t0, x0[0], dt);
        Real sigma2 = yProcess_->stdDeviation(t0, x0[1], dt);
        Real expa = std::exp(-a_*dt), expb = std::exp(-b_*dt);
        Real H = (rho_*sigma_*eta_)/(a_+b_)*(1-expa*expb);
        Real den =
            (0.5*sigma_*eta_)*std::sqrt((1-expa*expa)*(1-expb*expb)/(a_*b_));
        Real newRho = H/den;
        tmp[0][0] = sigma1;
        tmp[0][1] = 0.0;
        tmp[1][0] = newRho*sigma2;
        tmp[1][1] = std::sqrt(1.0-newRho*newRho)*sigma2;
        return tmp;
    }

    Matrix G2ForwardProcess::covariance(Time t0, const Array& x0, Time dt) const {
        Matrix sigma = stdDeviation(t0, x0, dt);
        Matrix result = sigma*transpose(sigma);
        return result;
    }

    Real G2ForwardProcess::xForwardDrift(Time t, Time T) const {
        Real expatT = std::exp(-a_*(T-t));
        Real expbtT = std::exp(-b_*(T-t));

        return -(sigma_*sigma_/a_) * (1-expatT)
              - (rho_*sigma_*eta_/b_) * (1-expbtT);
    }

    Real G2ForwardProcess::yForwardDrift(Time t, Time T) const {
        Real expatT = std::exp(-a_*(T-t));
        Real expbtT = std::exp(-b_*(T-t));

        return -(eta_*eta_/b_) * (1-expbtT)
              - (rho_*sigma_*eta_/a_) * (1-expatT);
    }

    Real G2ForwardProcess::Mx_T(Real s, Real t, Real T) const {
        Real M;
        M = ( (sigma_*sigma_)/(a_*a_) + (rho_*sigma_*eta_)/(a_*b_) )
          * (1-std::exp(-a_*(t-s)));
        M += -(sigma_*sigma_)/(2*a_*a_) *
              (std::exp(-a_*(T-t))-std::exp(-a_*(T+t-2*s)));
        M += -(rho_*sigma_*eta_)/(b_*(a_+b_))
            * (std::exp(-b_*(T-t)) -std::exp(-b_*T-a_*t+(a_+b_)*s));
        return M;
    }

    Real G2ForwardProcess::My_T(Real s, Real t, Real T) const {
        Real M;
        M = ( (eta_*eta_)/(b_*b_) + (rho_*sigma_*eta_)/(a_*b_) )
          * (1-std::exp(-b_*(t-s)));
        M += -(eta_*eta_)/(2*b_*b_) *
              (std::exp(-b_*(T-t))-std::exp(-b_*(T+t-2*s)));
        M += -(rho_*sigma_*eta_)/(a_*(a_+b_))
            * (std::exp(-a_*(T-t))-std::exp(-a_*T-b_*t+(a_+b_)*s));
        return M;
    }

}


Coverage Report

Created: 2025-08-28 06:30

Line	Count	Source (jump to first uncovered line)
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
23		namespace QuantLib {
24
25		G2Process::G2Process(Real a, Real sigma, Real b, Real eta, Real rho)
26	0	: a_(a), sigma_(sigma), b_(b), eta_(eta), rho_(rho),
27	0	xProcess_(new QuantLib::OrnsteinUhlenbeckProcess(a, sigma, 0.0)),
28	0	yProcess_(new QuantLib::OrnsteinUhlenbeckProcess(b, eta, 0.0)) {}
29
30	0	Size G2Process::size() const {
31	0	return 2;
32	0	}
33
34	0	Array G2Process::initialValues() const {
35	0	return { x0_, y0_ };
36	0	}
37
38	0	Array G2Process::drift(Time t, const Array& x) const {
39	0	return {
40	0	xProcess_->drift(t, x[0]),
41	0	yProcess_->drift(t, x[1])
42	0	};
43	0	}
44
45	0	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	0	Matrix tmp(2,2);
54	0	Real sigma1 = sigma_;
55	0	Real sigma2 = eta_;
56	0	tmp[0][0] = sigma1; tmp[0][1] = 0.0;
57	0	tmp[1][0] = rho_sigma1; tmp[1][1] = std::sqrt(1.0-rho_rho_)*sigma2;
58	0	return tmp;
59	0	}
60
61		Array G2Process::expectation(Time t0, const Array& x0,
62	0	Time dt) const {
63	0	return {
64	0	xProcess_->expectation(t0, x0[0], dt),
65	0	yProcess_->expectation(t0, x0[1], dt)
66	0	};
67	0	}
68
69	0	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	0	Matrix tmp(2,2);
78	0	Real sigma1 = xProcess_->stdDeviation(t0, x0[0], dt);
79	0	Real sigma2 = yProcess_->stdDeviation(t0, x0[1], dt);
80	0	Real expa = std::exp(-a_dt), expb = std::exp(-b_dt);
81	0	Real H = (rho_sigma_eta_)/(a_+b_)(1-expaexpb);
82	0	Real den =
83	0	(0.5sigma_eta_)std::sqrt((1-expaexpa)(1-expbexpb)/(a_*b_));
84	0	Real newRho = H/den;
85	0	tmp[0][0] = sigma1;
86	0	tmp[0][1] = 0.0;
87	0	tmp[1][0] = newRho*sigma2;
88	0	tmp[1][1] = std::sqrt(1.0-newRhonewRho)sigma2;
89	0	return tmp;
90	0	}
91
92	0	Matrix G2Process::covariance(Time t0, const Array& x0, Time dt) const {
93	0	Matrix sigma = stdDeviation(t0, x0, dt);
94	0	Matrix result = sigma*transpose(sigma);
95	0	return result;
96	0	}
97
98	0	Real G2Process::x0() const {
99	0	return x0_;
100	0	}
101
102	0	Real G2Process::y0() const {
103	0	return y0_;
104	0	}
105
106	0	Real G2Process::a() const {
107	0	return a_;
108	0	}
109
110	0	Real G2Process::sigma() const {
111	0	return sigma_;
112	0	}
113
114	0	Real G2Process::b() const {
115	0	return b_;
116	0	}
117
118	0	Real G2Process::eta() const {
119	0	return eta_;
120	0	}
121
122	0	Real G2Process::rho() const {
123	0	return rho_;
124	0	}
125
126
127		G2ForwardProcess::G2ForwardProcess(Real a, Real sigma, Real b, Real eta, Real rho)
128	0	: a_(a), sigma_(sigma), b_(b), eta_(eta), rho_(rho),
129	0	xProcess_(new QuantLib::OrnsteinUhlenbeckProcess(a, sigma, 0.0)),
130	0	yProcess_(new QuantLib::OrnsteinUhlenbeckProcess(b, eta, 0.0)) {}
131
132	0	Size G2ForwardProcess::size() const {
133	0	return 2;
134	0	}
135
136	0	Array G2ForwardProcess::initialValues() const {
137	0	return { x0_, y0_ };
138	0	}
139
140	0	Array G2ForwardProcess::drift(Time t, const Array& x) const {
141	0	return {
142	0	xProcess_->drift(t, x[0]) + xForwardDrift(t, T_),
143	0	yProcess_->drift(t, x[1]) + yForwardDrift(t, T_)
144	0	};
145	0	}
146
147	0	Matrix G2ForwardProcess::diffusion(Time, const Array&) const {
148	0	Matrix tmp(2,2);
149	0	Real sigma1 = sigma_;
150	0	Real sigma2 = eta_;
151	0	tmp[0][0] = sigma1; tmp[0][1] = 0.0;
152	0	tmp[1][0] = rho_sigma1; tmp[1][1] = std::sqrt(1.0-rho_rho_)*sigma2;
153	0	return tmp;
154	0	}
155
156		Array G2ForwardProcess::expectation(Time t0, const Array& x0,
157	0	Time dt) const {
158	0	return {
159	0	xProcess_->expectation(t0, x0[0], dt) - Mx_T(t0, t0+dt, T_),
160	0	yProcess_->expectation(t0, x0[1], dt) - My_T(t0, t0+dt, T_)
161	0	};
162	0	}
163
164	0	Matrix G2ForwardProcess::stdDeviation(Time t0, const Array& x0, Time dt) const {
165	0	Matrix tmp(2,2);
166	0	Real sigma1 = xProcess_->stdDeviation(t0, x0[0], dt);
167	0	Real sigma2 = yProcess_->stdDeviation(t0, x0[1], dt);
168	0	Real expa = std::exp(-a_dt), expb = std::exp(-b_dt);
169	0	Real H = (rho_sigma_eta_)/(a_+b_)(1-expaexpb);
170	0	Real den =
171	0	(0.5sigma_eta_)std::sqrt((1-expaexpa)(1-expbexpb)/(a_*b_));
172	0	Real newRho = H/den;
173	0	tmp[0][0] = sigma1;
174	0	tmp[0][1] = 0.0;
175	0	tmp[1][0] = newRho*sigma2;
176	0	tmp[1][1] = std::sqrt(1.0-newRhonewRho)sigma2;
177	0	return tmp;
178	0	}
179
180	0	Matrix G2ForwardProcess::covariance(Time t0, const Array& x0, Time dt) const {
181	0	Matrix sigma = stdDeviation(t0, x0, dt);
182	0	Matrix result = sigma*transpose(sigma);
183	0	return result;
184	0	}
185
186	0	Real G2ForwardProcess::xForwardDrift(Time t, Time T) const {
187	0	Real expatT = std::exp(-a_*(T-t));
188	0	Real expbtT = std::exp(-b_*(T-t));
189
190	0	return -(sigma_sigma_/a_) (1-expatT)
191	0	- (rho_sigma_eta_/b_) * (1-expbtT);
192	0	}
193
194	0	Real G2ForwardProcess::yForwardDrift(Time t, Time T) const {
195	0	Real expatT = std::exp(-a_*(T-t));
196	0	Real expbtT = std::exp(-b_*(T-t));
197
198	0	return -(eta_eta_/b_) (1-expbtT)
199	0	- (rho_sigma_eta_/a_) * (1-expatT);
200	0	}
201
202	0	Real G2ForwardProcess::Mx_T(Real s, Real t, Real T) const {
203	0	Real M;
204	0	M = ( (sigma_sigma_)/(a_a_) + (rho_sigma_eta_)/(a_*b_) )
205	0	* (1-std::exp(-a_*(t-s)));
206	0	M += -(sigma_sigma_)/(2a_a_)
207	0	(std::exp(-a_(T-t))-std::exp(-a_(T+t-2*s)));
208	0	M += -(rho_sigma_eta_)/(b_*(a_+b_))
209	0	* (std::exp(-b_(T-t)) -std::exp(-b_T-a_t+(a_+b_)s));
210	0	return M;
211	0	}
212
213	0	Real G2ForwardProcess::My_T(Real s, Real t, Real T) const {
214	0	Real M;
215	0	M = ( (eta_eta_)/(b_b_) + (rho_sigma_eta_)/(a_*b_) )
216	0	* (1-std::exp(-b_*(t-s)));
217	0	M += -(eta_eta_)/(2b_b_)
218	0	(std::exp(-b_(T-t))-std::exp(-b_(T+t-2*s)));
219	0	M += -(rho_sigma_eta_)/(a_*(a_+b_))
220	0	* (std::exp(-a_(T-t))-std::exp(-a_T-b_t+(a_+b_)s));
221	0	return M;
222	0	}
223
224		}
225