8478223765

Committed 29 Mar 2024 06:56AM UTC coverage: 72.492% (-0.005%) from 72.497%

Build # 8478223765

Build Type

Pull #1937

github

Committed by

web-flow

Commit Message

Merge 7a21fd833 into 9ff05427c

Pull Request Pull Request #1937: Added Interest-rate term structure in g2process constructor

Run Details

0 of 11 new or added lines in 1 file covered. (0.0%)

2 existing lines in 1 file now uncovered.

54966 of 75824 relevant lines covered (72.49%)

8679241.9 hits per line

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

0.0

/ql/processes/g2process.cpp

/* -*- 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(const Handle<YieldTermStructure>& h, 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)), h_(h) {}

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

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

    Array G2Process::drift(Time t, const Array& x) const {
        Real alpha_drift = sigma_*sigma_/(2*a_)*(1-std::exp(-2*a_*t));
        Real shift = 0.0001;
        Real f = h_->forwardRate(t, t, Continuous, NoFrequency);
        Real fup = h_->forwardRate(t+shift, t+shift, Continuous, NoFrequency);
        Real f_prime = (fup-f)/shift;
        alpha_drift += a_*f+f_prime;

        return {
            xProcess_->drift(t, x[0]) + alpha_drift,
            yProcess_->drift(t, x[1]) + alpha_drift
        };
    }

    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(const Handle<YieldTermStructure>& h, 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)),h_(h) {}

    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;
    }

}


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
NEW 25	G2Process::G2Process(const Handle<YieldTermStructure>& h, 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)),	×
NEW 28	yProcess_(new QuantLib::OrnsteinUhlenbeckProcess(b, eta, 0.0)), h_(h) {}	×
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 {	×
NEW 39	Real alpha_drift = sigma_sigma_/(2a_)(1-std::exp(-2a_*t));	×
40	Real shift = 0.0001;
NEW 41	Real f = h_->forwardRate(t, t, Continuous, NoFrequency);	×
NEW 42	Real fup = h_->forwardRate(t+shift, t+shift, Continuous, NoFrequency);	×
NEW 43	Real f_prime = (fup-f)/shift;	×
NEW 44	alpha_drift += a_*f+f_prime;	×
45
46	return {
NEW 47	xProcess_->drift(t, x[0]) + alpha_drift,	×
NEW 48	yProcess_->drift(t, x[1]) + alpha_drift	×
UNCOV 49	};	×
50	}
51
52	Matrix G2Process::diffusion(Time, const Array&) const {	×
53	/* the correlation matrix is
54	\| 1 rho \|
55	\| rho 1 \|
56	whose square root (which is used here) is
57	\| 1 0 \|
58	\| rho sqrt(1-rho^2) \|
59	*/
60	Matrix tmp(2,2);
61	Real sigma1 = sigma_;	×
62	Real sigma2 = eta_;	×
63	tmp[0][0] = sigma1; tmp[0][1] = 0.0;	×
64	tmp[1][0] = rho_sigma1; tmp[1][1] = std::sqrt(1.0-rho_rho_)*sigma2;	×
65	return tmp;	×
66	}
67
68	Array G2Process::expectation(Time t0, const Array& x0,	×
69	Time dt) const {
70	return {
71	xProcess_->expectation(t0, x0[0], dt),	×
72	yProcess_->expectation(t0, x0[1], dt)	×
73	};	×
74	}
75
76	Matrix G2Process::stdDeviation(Time t0, const Array& x0, Time dt) const {	×
77	/* the correlation matrix is
78	\| 1 rho \|
79	\| rho 1 \|
80	whose square root (which is used here) is
81	\| 1 0 \|
82	\| rho sqrt(1-rho^2) \|
83	*/
84	Matrix tmp(2,2);
85	Real sigma1 = xProcess_->stdDeviation(t0, x0[0], dt);	×
86	Real sigma2 = yProcess_->stdDeviation(t0, x0[1], dt);	×
87	Real expa = std::exp(-a_dt), expb = std::exp(-b_dt);	×
88	Real H = (rho_sigma_eta_)/(a_+b_)(1-expaexpb);	×
89	Real den =
90	(0.5sigma_eta_)std::sqrt((1-expaexpa)(1-expbexpb)/(a_*b_));	×
91	Real newRho = H/den;	×
92	tmp[0][0] = sigma1;	×
93	tmp[0][1] = 0.0;	×
94	tmp[1][0] = newRho*sigma2;	×
95	tmp[1][1] = std::sqrt(1.0-newRhonewRho)sigma2;	×
96	return tmp;	×
97	}
98
99	Matrix G2Process::covariance(Time t0, const Array& x0, Time dt) const {	×
100	Matrix sigma = stdDeviation(t0, x0, dt);	×
101	Matrix result = sigma*transpose(sigma);	×
102	return result;	×
103	}
104
105	Real G2Process::x0() const {	×
106	return x0_;	×
107	}
108
109	Real G2Process::y0() const {	×
110	return y0_;	×
111	}
112
113	Real G2Process::a() const {	×
114	return a_;	×
115	}
116
117	Real G2Process::sigma() const {	×
118	return sigma_;	×
119	}
120
121	Real G2Process::b() const {	×
122	return b_;	×
123	}
124
125	Real G2Process::eta() const {	×
126	return eta_;	×
127	}
128
129	Real G2Process::rho() const {	×
130	return rho_;	×
131	}
132
133
NEW 134	G2ForwardProcess::G2ForwardProcess(const Handle<YieldTermStructure>& h, Real a, Real sigma, Real b, Real eta, Real rho)	×
135	: a_(a), sigma_(sigma), b_(b), eta_(eta), rho_(rho),	×
136	xProcess_(new QuantLib::OrnsteinUhlenbeckProcess(a, sigma, 0.0)),	×
NEW 137	yProcess_(new QuantLib::OrnsteinUhlenbeckProcess(b, eta, 0.0)),h_(h) {}	×
138
139	Size G2ForwardProcess::size() const {	×
140	return 2;	×
141	}
142
143	Array G2ForwardProcess::initialValues() const {	×
144	return { x0_, y0_ };	×
145	}
146
147	Array G2ForwardProcess::drift(Time t, const Array& x) const {	×
148	return {
149	xProcess_->drift(t, x[0]) + xForwardDrift(t, T_),	×
150	yProcess_->drift(t, x[1]) + yForwardDrift(t, T_)	×
151	};	×
152	}
153
154	Matrix G2ForwardProcess::diffusion(Time, const Array&) const {	×
155	Matrix tmp(2,2);
156	Real sigma1 = sigma_;	×
157	Real sigma2 = eta_;	×
158	tmp[0][0] = sigma1; tmp[0][1] = 0.0;	×
159	tmp[1][0] = rho_sigma1; tmp[1][1] = std::sqrt(1.0-rho_rho_)*sigma2;	×
160	return tmp;	×
161	}
162
163	Array G2ForwardProcess::expectation(Time t0, const Array& x0,	×
164	Time dt) const {
165	return {
166	xProcess_->expectation(t0, x0[0], dt) - Mx_T(t0, t0+dt, T_),	×
167	yProcess_->expectation(t0, x0[1], dt) - My_T(t0, t0+dt, T_)	×
168	};	×
169	}
170
171	Matrix G2ForwardProcess::stdDeviation(Time t0, const Array& x0, Time dt) const {	×
172	Matrix tmp(2,2);
173	Real sigma1 = xProcess_->stdDeviation(t0, x0[0], dt);	×
174	Real sigma2 = yProcess_->stdDeviation(t0, x0[1], dt);	×
175	Real expa = std::exp(-a_dt), expb = std::exp(-b_dt);	×
176	Real H = (rho_sigma_eta_)/(a_+b_)(1-expaexpb);	×
177	Real den =
178	(0.5sigma_eta_)std::sqrt((1-expaexpa)(1-expbexpb)/(a_*b_));	×
179	Real newRho = H/den;	×
UNCOV 180	tmp[0][0] = sigma1;	×
181	tmp[0][1] = 0.0;	×
182	tmp[1][0] = newRho*sigma2;	×
183	tmp[1][1] = std::sqrt(1.0-newRhonewRho)sigma2;	×
184	return tmp;	×
185	}
186
187	Matrix G2ForwardProcess::covariance(Time t0, const Array& x0, Time dt) const {	×
188	Matrix sigma = stdDeviation(t0, x0, dt);	×
189	Matrix result = sigma*transpose(sigma);	×
190	return result;	×
191	}
192
193	Real G2ForwardProcess::xForwardDrift(Time t, Time T) const {	×
194	Real expatT = std::exp(-a_*(T-t));	×
195	Real expbtT = std::exp(-b_*(T-t));	×
196
197	return -(sigma_sigma_/a_) (1-expatT)	×
198	- (rho_sigma_eta_/b_) * (1-expbtT);	×
199	}
200
201	Real G2ForwardProcess::yForwardDrift(Time t, Time T) const {	×
202	Real expatT = std::exp(-a_*(T-t));	×
203	Real expbtT = std::exp(-b_*(T-t));	×
204
205	return -(eta_eta_/b_) (1-expbtT)	×
206	- (rho_sigma_eta_/a_) * (1-expatT);	×
207	}
208
209	Real G2ForwardProcess::Mx_T(Real s, Real t, Real T) const {	×
210	Real M;
211	M = ( (sigma_sigma_)/(a_a_) + (rho_sigma_eta_)/(a_*b_) )	×
212	* (1-std::exp(-a_*(t-s)));	×
213	M += -(sigma_sigma_)/(2a_a_)	×
214	(std::exp(-a_(T-t))-std::exp(-a_(T+t-2*s)));	×
215	M += -(rho_sigma_eta_)/(b_*(a_+b_))	×
216	* (std::exp(-b_(T-t)) -std::exp(-b_T-a_t+(a_+b_)s));	×
217	return M;	×
218	}
219
220	Real G2ForwardProcess::My_T(Real s, Real t, Real T) const {	×
221	Real M;
222	M = ( (eta_eta_)/(b_b_) + (rho_sigma_eta_)/(a_*b_) )	×
223	* (1-std::exp(-b_*(t-s)));	×
224	M += -(eta_eta_)/(2b_b_)	×
225	(std::exp(-b_(T-t))-std::exp(-b_(T+t-2*s)));	×
226	M += -(rho_sigma_eta_)/(a_*(a_+b_))	×
227	* (std::exp(-a_(T-t))-std::exp(-a_T-b_t+(a_+b_)s));	×
228	return M;	×
229	}
230
231	}
232

lballabio / QuantLib / 8478223765

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