QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
gaussian1dmodel.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) 2013, 2015 Peter Caspers
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
22#include <ql/payoff.hpp>
23#include <cmath>
24
25namespace QuantLib {
26
28 const Date& referenceDate,
29 const Real y,
30 const ext::shared_ptr<IborIndex>& iborIdx) const {
31
32 QL_REQUIRE(iborIdx != nullptr, "no ibor index given");
33
34 calculate();
35
36 if (fixing <= (evaluationDate_ + (enforcesTodaysHistoricFixings_ ? 0 : -1)))
37 return iborIdx->fixing(fixing);
38
39 Handle<YieldTermStructure> yts = iborIdx->forwardingTermStructure(); // might be empty, then
40 // use model curve
41
42 Date valueDate = iborIdx->valueDate(fixing);
43 Date endDate = iborIdx->fixingCalendar().advance(
44 valueDate, iborIdx->tenor(), iborIdx->businessDayConvention(), iborIdx->endOfMonth());
45 // FIXME Here we should use the calculation date calendar ?
46 Real dcf = iborIdx->dayCounter().yearFraction(valueDate, endDate);
47
48 return (zerobond(valueDate, referenceDate, y, yts) -
49 zerobond(endDate, referenceDate, y, yts)) /
50 (dcf * zerobond(endDate, referenceDate, y, yts));
51}
52
54 const Period& tenor,
55 const Date& referenceDate,
56 const Real y,
57 const ext::shared_ptr<SwapIndex>& swapIdx) const {
58
59 QL_REQUIRE(swapIdx != nullptr, "no swap index given");
60
61 calculate();
62
63 if (fixing <= (evaluationDate_ + (enforcesTodaysHistoricFixings_ ? 0 : -1)))
64 return swapIdx->fixing(fixing);
65
67 swapIdx->iborIndex()->forwardingTermStructure();
69 swapIdx->discountingTermStructure(); // either might be empty, then
70 // use model curve
71
72 Schedule sched, floatSched;
73
74 ext::shared_ptr<VanillaSwap> underlying =
75 underlyingSwap(swapIdx, fixing, tenor);
76
77 sched = underlying->fixedSchedule();
78
79 ext::shared_ptr<OvernightIndexedSwapIndex> oisIdx =
80 ext::dynamic_pointer_cast<OvernightIndexedSwapIndex>(swapIdx);
81 if (oisIdx != nullptr) {
82 floatSched = sched;
83 } else {
84 floatSched = underlying->floatingSchedule();
85 }
86
87 Real annuity = swapAnnuity(fixing, tenor, referenceDate, y,
88 swapIdx); // should be fine for
89 // overnightindexed swap indices as
90 // well
91 Rate floatleg = 0.0;
92 if (ytsf.empty() && ytsd.empty()) { // simple 100-formula can be used
93 // only in one curve setup
94 floatleg =
95 (zerobond(sched.dates().front(), referenceDate, y) -
96 zerobond(sched.calendar().adjust(sched.dates().back(),
97 underlying->paymentConvention()),
98 referenceDate, y));
99 } else {
100 for (Size i = 1; i < floatSched.size(); i++) {
101 floatleg +=
102 (zerobond(floatSched[i - 1], referenceDate, y, ytsf) /
103 zerobond(floatSched[i], referenceDate, y, ytsf) -
104 1.0) *
105 zerobond(floatSched.calendar().adjust(
106 floatSched[i], underlying->paymentConvention()),
107 referenceDate, y, ytsd);
108 }
109 }
110 return floatleg / annuity;
111}
112
114 const Period& tenor,
115 const Date& referenceDate,
116 const Real y,
117 const ext::shared_ptr<SwapIndex>& swapIdx) const {
118
119 QL_REQUIRE(swapIdx != nullptr, "no swap index given");
120
121 calculate();
122
124 swapIdx->discountingTermStructure(); // might be empty, then use
125 // model curve
126
127 ext::shared_ptr<VanillaSwap> underlying =
128 underlyingSwap(swapIdx, fixing, tenor);
129
130 Schedule sched = underlying->fixedSchedule();
131
132 Real annuity = 0.0;
133 for (unsigned int j = 1; j < sched.size(); j++) {
134 annuity += zerobond(sched.calendar().adjust(
135 sched.date(j), underlying->paymentConvention()),
136 referenceDate, y, ytsd) *
137 swapIdx->dayCounter().yearFraction(sched.date(j - 1),
138 sched.date(j));
139 }
140 return annuity;
141}
142
144 const Option::Type &type, const Date &expiry, const Date &valueDate,
145 const Date &maturity, const Rate strike, const Date &referenceDate,
146 const Real y, const Handle<YieldTermStructure> &yts, const Real yStdDevs,
147 const Size yGridPoints, const bool extrapolatePayoff,
148 const bool flatPayoffExtrapolation) const {
149
150 calculate();
151
152 Time fixingTime = termStructure()->timeFromReference(expiry);
153 Time referenceTime =
154 referenceDate == Null<Date>()
155 ? 0.0
156 : termStructure()->timeFromReference(referenceDate);
157
158 Array yg = yGrid(yStdDevs, yGridPoints, fixingTime, referenceTime, y);
159 Array z = yGrid(yStdDevs, yGridPoints);
160
161 Array p(yg.size());
162
163 for (Size i = 0; i < yg.size(); i++) {
164 Real expValDsc = zerobond(valueDate, expiry, yg[i], yts);
165 Real discount =
166 zerobond(maturity, expiry, yg[i], yts) / expValDsc;
167 p[i] =
168 std::max((type == Option::Call ? 1.0 : -1.0) * (discount - strike),
169 0.0) /
170 numeraire(fixingTime, yg[i], yts) * expValDsc;
171 }
172
174 z.begin(), z.end(), p.begin(), CubicInterpolation::Spline, true,
176
177 Real price = 0.0;
178 for (Size i = 0; i < z.size() - 1; i++) {
180 0.0, payoff.cCoefficients()[i], payoff.bCoefficients()[i],
181 payoff.aCoefficients()[i], p[i], z[i], z[i], z[i + 1]);
182 }
183 if (extrapolatePayoff) {
184 if (flatPayoffExtrapolation) {
186 0.0, 0.0, 0.0, 0.0, p[z.size() - 2], z[z.size() - 2],
187 z[z.size() - 1], 100.0);
188 price += gaussianShiftedPolynomialIntegral(0.0, 0.0, 0.0, 0.0, p[0],
189 z[0], -100.0, z[0]);
190 } else {
191 if (type == Option::Call)
193 0.0, payoff.cCoefficients()[z.size() - 2],
194 payoff.bCoefficients()[z.size() - 2],
195 payoff.aCoefficients()[z.size() - 2], p[z.size() - 2],
196 z[z.size() - 2], z[z.size() - 1], 100.0);
197 if (type == Option::Put)
199 0.0, payoff.cCoefficients()[0], payoff.bCoefficients()[0],
200 payoff.aCoefficients()[0], p[0], z[0], -100.0, z[0]);
201 }
202 }
203
204 return numeraire(referenceTime, y, yts) * price;
205}
206
208 const Real a, const Real b, const Real c, const Real d, const Real e,
209 const Real y0, const Real y1) {
210
211#ifdef GAUSS1D_ENABLE_NTL
212 const boost::math::ntl::RR aa = 4.0 * a, ba = 2.0 * M_SQRT2 * b,
213 ca = 2.0 * c, da = M_SQRT2 * d;
214 const boost::math::ntl::RR x0 = y0 * M_SQRT1_2, x1 = y1 * M_SQRT1_2;
215 const boost::math::ntl::RR res =
216 (0.125 * (3.0 * aa + 2.0 * ca + 4.0 * e) * std::erf(x1) -
217 1.0 / (4.0 * M_SQRTPI) * std::exp(-x1 * x1) *
218 (2.0 * aa * x1 * x1 * x1 + 3.0 * aa * x1 +
219 2.0 * ba * (x1 * x1 + 1.0) + 2.0 * ca * x1 + 2.0 * da)) -
220 (0.125 * (3.0 * aa + 2.0 * ca + 4.0 * e) * std::erf(x0) -
221 1.0 / (4.0 * M_SQRTPI) * std::exp(-x0 * x0) *
222 (2.0 * aa * x0 * x0 * x0 + 3.0 * aa * x0 +
223 2.0 * ba * (x0 * x0 + 1.0) + 2.0 * ca * x0 + 2.0 * da));
224 return NTL::to_double(res.value());
225#else
226 const Real aa = 4.0 * a, ba = 2.0 * M_SQRT2 * b, ca = 2.0 * c,
227 da = M_SQRT2 * d;
228 const Real x0 = y0 * M_SQRT1_2, x1 = y1 * M_SQRT1_2;
229 return (0.125 * (3.0 * aa + 2.0 * ca + 4.0 * e) * std::erf(x1) -
230 1.0 / (4.0 * M_SQRTPI) * std::exp(-x1 * x1) *
231 (2.0 * aa * x1 * x1 * x1 + 3.0 * aa * x1 +
232 2.0 * ba * (x1 * x1 + 1.0) + 2.0 * ca * x1 + 2.0 * da)) -
233 (0.125 * (3.0 * aa + 2.0 * ca + 4.0 * e) * std::erf(x0) -
234 1.0 / (4.0 * M_SQRTPI) * std::exp(-x0 * x0) *
235 (2.0 * aa * x0 * x0 * x0 + 3.0 * aa * x0 +
236 2.0 * ba * (x0 * x0 + 1.0) + 2.0 * ca * x0 + 2.0 * da));
237#endif
238}
239
241 const Real a, const Real b, const Real c, const Real d, const Real e,
242 const Real h, const Real x0, const Real x1) {
244 a, -4.0 * a * h + b, 6.0 * a * h * h - 3.0 * b * h + c,
245 -4 * a * h * h * h + 3.0 * b * h * h - 2.0 * c * h + d,
246 a * h * h * h * h - b * h * h * h + c * h * h - d * h + e, x0, x1);
247}
248
250 const Real stdDevs, const int gridPoints, const Real T, const Real t, const Real y) const {
251
252 // we use that the standard deviation is independent of $x$ here !
253
254 QL_REQUIRE(stateProcess_ != nullptr, "state process not set");
255
256 Array result(2 * gridPoints + 1, 0.0);
257
258 Real x_t, e_0_t, e_t_T, stdDev_0_t, stdDev_t_T;
259 Real stdDev_0_T = stateProcess_->stdDeviation(0.0, 0.0, T);
260 Real e_0_T = stateProcess_->expectation(0.0, 0.0, T);
261
262 if (t < QL_EPSILON) {
263 // stdDev_0_t = 0.0;
264 stdDev_t_T = stdDev_0_T;
265 // e_0_t = 0.0;
266 // x_t = 0.0;
267 e_t_T = e_0_T;
268 } else {
269 stdDev_0_t = stateProcess_->stdDeviation(0.0, 0.0, t);
270 stdDev_t_T = stateProcess_->stdDeviation(t, 0.0, T - t);
271 e_0_t = stateProcess_->expectation(0.0, 0.0, t);
272 x_t = y * stdDev_0_t + e_0_t;
273 e_t_T = stateProcess_->expectation(t, x_t, T - t);
274 }
275
276 Real h = stdDevs / ((Real)gridPoints);
277
278 for (int j = -gridPoints; j <= gridPoints; j++) {
279 result[j + gridPoints] =
280 (e_t_T + stdDev_t_T * ((Real)j) * h - e_0_T) / stdDev_0_T;
281 }
282
283 return result;
284}
285}
1-D array used in linear algebra.
Definition: array.hpp:52
const_iterator end() const
Definition: array.hpp:511
Size size() const
dimension of the array
Definition: array.hpp:495
const_iterator begin() const
Definition: array.hpp:503
Date adjust(const Date &, BusinessDayConvention convention=Following) const
Definition: calendar.cpp:84
Cubic interpolation between discrete points.
Concrete date class.
Definition: date.hpp:125
static Date endOfMonth(const Date &d)
last day of the month to which the given date belongs
Definition: date.hpp:428
static Date advance(const Date &d, Integer units, TimeUnit)
Definition: date.cpp:139
static Real gaussianPolynomialIntegral(Real a, Real b, Real c, Real d, Real e, Real x0, Real x1)
static Real gaussianShiftedPolynomialIntegral(Real a, Real b, Real c, Real d, Real e, Real h, Real x0, Real x1)
ext::shared_ptr< VanillaSwap > underlyingSwap(const ext::shared_ptr< SwapIndex > &index, const Date &expiry, const Period &tenor) const
Real numeraire(Time t, Real y=0.0, const Handle< YieldTermStructure > &yts=Handle< YieldTermStructure >()) const
Real swapRate(const Date &fixing, const Period &tenor, const Date &referenceDate=Null< Date >(), Real y=0.0, const ext::shared_ptr< SwapIndex > &swapIdx=ext::shared_ptr< SwapIndex >()) const
Array yGrid(Real yStdDevs, int gridPoints, Real T=1.0, Real t=0, Real y=0) const
Real forwardRate(const Date &fixing, const Date &referenceDate=Null< Date >(), Real y=0.0, const ext::shared_ptr< IborIndex > &iborIdx=ext::shared_ptr< IborIndex >()) const
Real zerobond(Time T, Time t=0.0, Real y=0.0, const Handle< YieldTermStructure > &yts=Handle< YieldTermStructure >()) const
Real zerobondOption(const Option::Type &type, const Date &expiry, const Date &valueDate, const Date &maturity, Rate strike, const Date &referenceDate=Null< Date >(), Real y=0.0, const Handle< YieldTermStructure > &yts=Handle< YieldTermStructure >(), Real yStdDevs=7.0, Size yGridPoints=64, bool extrapolatePayoff=true, bool flatPayoffExtrapolation=false) const
Real swapAnnuity(const Date &fixing, const Period &tenor, const Date &referenceDate=Null< Date >(), Real y=0.0, const ext::shared_ptr< SwapIndex > &swapIdx=ext::shared_ptr< SwapIndex >()) const
ext::shared_ptr< StochasticProcess1D > stateProcess_
Shared handle to an observable.
Definition: handle.hpp:41
bool empty() const
checks if the contained shared pointer points to anything
Definition: handle.hpp:191
virtual void calculate() const
Definition: lazyobject.hpp:253
template class providing a null value for a given type.
Definition: null.hpp:76
Payment schedule.
Definition: schedule.hpp:40
const Calendar & calendar() const
Definition: schedule.hpp:176
const std::vector< Date > & dates() const
Definition: schedule.hpp:75
const Date & date(Size i) const
Definition: schedule.hpp:160
Size size() const
Definition: schedule.hpp:69
const Handle< YieldTermStructure > & termStructure() const
Definition: model.hpp:77
cubic interpolation between discrete points
const DefaultType & t
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
Definition: errors.hpp:117
Date d
ext::function< Real(Real)> b
basic interface for one factor interest rate models
#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 Rate
interest rates
Definition: types.hpp:70
std::size_t Size
size of a container
Definition: types.hpp:58
ext::shared_ptr< QuantLib::Payoff > payoff
#define M_SQRTPI
#define M_SQRT2
#define M_SQRT1_2
Definition: any.hpp:35
Option payoff classes.