QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
xabrinterpolation.hpp
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) 2006 Ferdinando Ametrano
5 Copyright (C) 2007 Marco Bianchetti
6 Copyright (C) 2007 François du Vignaud
7 Copyright (C) 2007 Giorgio Facchinetti
8 Copyright (C) 2006 Mario Pucci
9 Copyright (C) 2006 StatPro Italia srl
10 Copyright (C) 2014 Peter Caspers
11
12 This file is part of QuantLib, a free-software/open-source library
13 for financial quantitative analysts and developers - http://quantlib.org/
14
15 QuantLib is free software: you can redistribute it and/or modify it
16 under the terms of the QuantLib license. You should have received a
17 copy of the license along with this program; if not, please email
18 <quantlib-dev@lists.sf.net>. The license is also available online at
19 <http://quantlib.org/license.shtml>.
20
21 This program is distributed in the hope that it will be useful, but WITHOUT
22 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
23 FOR A PARTICULAR PURPOSE. See the license for more details.
24*/
25
26/*! \file xabrinterpolation.hpp
27 \brief generic interpolation class for sabr style underlying models
28 like the Hagan 2002 expansion, Doust's no arbitrage sabr,
29 Andreasen's zabr expansion for the masses and similar
30*/
31
32#ifndef ql_xabr_interpolation_hpp
33#define ql_xabr_interpolation_hpp
34
45#include <ql/utilities/null.hpp>
46#include <utility>
47
48namespace QuantLib {
49
50namespace detail {
51
52template <typename Model> class XABRCoeffHolder {
53 public:
55 const Real& forward,
56 const std::vector<Real>& params,
57 const std::vector<bool>& paramIsFixed,
58 std::vector<Real> addParams)
59 : t_(t), forward_(forward), params_(params), paramIsFixed_(paramIsFixed.size(), false),
60 error_(Null<Real>()), maxError_(Null<Real>()), addParams_(std::move(addParams)) {
61 QL_REQUIRE(t > 0.0, "expiry time must be positive: " << t
62 << " not allowed");
63 QL_REQUIRE(params.size() == Model().dimension(),
64 "wrong number of parameters (" << params.size()
65 << "), should be "
66 << Model().dimension());
67 QL_REQUIRE(paramIsFixed.size() == Model().dimension(),
68 "wrong number of fixed parameters flags ("
69 << paramIsFixed.size() << "), should be "
70 << Model().dimension());
71
72 for (Size i = 0; i < params.size(); ++i) {
73 if (params[i] != Null<Real>())
74 paramIsFixed_[i] = paramIsFixed[i];
75 }
76 Model().defaultValues(params_, paramIsFixed_, forward_, t_, addParams_);
78 }
79 virtual ~XABRCoeffHolder() = default;
80
82 modelInstance_ = Model().instance(t_, forward_, params_, addParams_);
83 }
84
85 /*! Expiry, Forward */
87 const Real &forward_;
88 /*! Parameters */
89 std::vector<Real> params_;
90 std::vector<bool> paramIsFixed_;
91 std::vector<Real> weights_;
92 /*! Interpolation results */
95 /*! Model instance (if required) */
96 ext::shared_ptr<typename Model::type> modelInstance_;
97 /*! additional parameters */
98 std::vector<Real> addParams_;
99};
100
101template <class I1, class I2, typename Model>
103 public XABRCoeffHolder<Model> {
104 public:
105 XABRInterpolationImpl(const I1& xBegin,
106 const I1& xEnd,
107 const I2& yBegin,
108 Time t,
109 const Real& forward,
110 const std::vector<Real>& params,
111 const std::vector<bool>& paramIsFixed,
112 bool vegaWeighted,
113 ext::shared_ptr<EndCriteria> endCriteria,
114 ext::shared_ptr<OptimizationMethod> optMethod,
115 const Real errorAccept,
116 const bool useMaxError,
117 const Size maxGuesses,
118 const std::vector<Real>& addParams = std::vector<Real>(),
120 : Interpolation::templateImpl<I1, I2>(xBegin, xEnd, yBegin, 1),
121 XABRCoeffHolder<Model>(t, forward, params, paramIsFixed, addParams),
122 endCriteria_(std::move(endCriteria)), optMethod_(std::move(optMethod)),
123 errorAccept_(errorAccept), useMaxError_(useMaxError), maxGuesses_(maxGuesses),
124 vegaWeighted_(vegaWeighted), volatilityType_(volatilityType) {
125 // if no optimization method or endCriteria is provided, we provide one
126 if (!optMethod_)
127 optMethod_ = ext::shared_ptr<OptimizationMethod>(
128 new LevenbergMarquardt(1e-8, 1e-8, 1e-8));
129 // optMethod_ = ext::shared_ptr<OptimizationMethod>(new
130 // Simplex(0.01));
131 if (!endCriteria_) {
132 endCriteria_ = ext::make_shared<EndCriteria>(
133 60000, 100, 1e-8, 1e-8, 1e-8);
134 }
135 this->weights_ =
136 std::vector<Real>(xEnd - xBegin, 1.0 / (xEnd - xBegin));
137 }
138
139 void update() override {
140
141 this->updateModelInstance();
142
143 // we should also check that y contains positive values only
144
145 // we must update weights if it is vegaWeighted
146 if (vegaWeighted_) {
147 I1 x = this->xBegin_;
148 I2 y = this->yBegin_;
149 // std::vector<Real>::iterator w = weights_.begin();
150 this->weights_.clear();
151 Real weightsSum = 0.0;
152 for (; x != this->xEnd_; ++x, ++y) {
153 Real stdDev = std::sqrt((*y) * (*y) * this->t_);
154 this->weights_.push_back(Model().weight(*x, this->forward_, stdDev,
155 this->addParams_));
156 weightsSum += this->weights_.back();
157 }
158 // weight normalization
159 auto w = this->weights_.begin();
160 for (; w != this->weights_.end(); ++w)
161 *w /= weightsSum;
162 }
163
164 // there is nothing to optimize
165 if (std::accumulate(this->paramIsFixed_.begin(),
166 this->paramIsFixed_.end(), true,
167 std::logical_and<>())) {
168 this->error_ = interpolationError();
171 return;
172 } else {
173 XABRError costFunction(this);
174
175 Array guess(Model().dimension());
176 for (Size i = 0; i < guess.size(); ++i)
177 guess[i] = this->params_[i];
178
179 Size iterations = 0;
180 Size freeParameters = 0;
181 Real bestError = QL_MAX_REAL;
182 Array bestParameters;
183 for (Size i = 0; i < Model().dimension(); ++i)
184 if (!this->paramIsFixed_[i])
185 ++freeParameters;
186 HaltonRsg halton(freeParameters, 42);
187 EndCriteria::Type tmpEndCriteria;
188 Real tmpInterpolationError;
189
190 do {
191
192 if (iterations > 0) {
193 const auto& s = halton.nextSequence();
194 Model().guess(guess, this->paramIsFixed_, this->forward_,
195 this->t_, s.value, this->addParams_);
196 for (Size i = 0; i < this->paramIsFixed_.size(); ++i)
197 if (this->paramIsFixed_[i])
198 guess[i] = this->params_[i];
199 }
200
201 Array inversedTransformatedGuess(Model().inverse(
202 guess, this->paramIsFixed_, this->params_, this->forward_));
203
204 ProjectedCostFunction constrainedXABRError(
205 costFunction, inversedTransformatedGuess,
206 this->paramIsFixed_);
207
208 Array projectedGuess(
209 constrainedXABRError.project(inversedTransformatedGuess));
210
211 NoConstraint constraint;
212 Problem problem(constrainedXABRError, constraint,
213 projectedGuess);
214 tmpEndCriteria = optMethod_->minimize(problem, *endCriteria_);
215 Array projectedResult(problem.currentValue());
216 Array transfResult(
217 constrainedXABRError.include(projectedResult));
218
219 Array result = Model().direct(transfResult, this->paramIsFixed_,
220 this->params_, this->forward_);
221 tmpInterpolationError = useMaxError_ ? interpolationMaxError()
223
224 if (tmpInterpolationError < bestError) {
225 bestError = tmpInterpolationError;
226 bestParameters = result;
227 this->XABREndCriteria_ = tmpEndCriteria;
228 }
229
230 } while (++iterations < maxGuesses_ &&
231 tmpInterpolationError > errorAccept_);
232
233 for (Size i = 0; i < bestParameters.size(); ++i)
234 this->params_[i] = bestParameters[i];
235
236 this->error_ = interpolationError();
238 }
239 }
240
241 Real value(Real x) const override { return this->modelInstance_->volatility(x, volatilityType_); }
242
243 Real primitive(Real) const override { QL_FAIL("XABR primitive not implemented"); }
244 Real derivative(Real) const override { QL_FAIL("XABR derivative not implemented"); }
245 Real secondDerivative(Real) const override { QL_FAIL("XABR secondDerivative not implemented"); }
246
247 // calculate total squared weighted difference (L2 norm)
249 Real error, totalError = 0.0;
250 I1 x = this->xBegin_;
251 I2 y = this->yBegin_;
252 auto w = this->weights_.begin();
253 for (; x != this->xEnd_; ++x, ++y, ++w) {
254 error = (value(*x) - *y);
255 totalError += error * error * (*w);
256 }
257 return totalError;
258 }
259
260 // calculate weighted differences
262 Array results(this->xEnd_ - this->xBegin_);
263 I1 x = this->xBegin_;
264 Array::iterator r = results.begin();
265 I2 y = this->yBegin_;
266 auto w = this->weights_.begin();
267 for (; x != this->xEnd_; ++x, ++r, ++w, ++y) {
268 *r = (value(*x) - *y) * std::sqrt(*w);
269 }
270 return results;
271 }
272
274 Size n = this->xEnd_ - this->xBegin_;
275 Real squaredError = interpolationSquaredError();
276 return std::sqrt(n * squaredError / (n==1 ? 1 : (n - 1)));
277 }
278
280 Real error, maxError = QL_MIN_REAL;
281 I1 i = this->xBegin_;
282 I2 j = this->yBegin_;
283 for (; i != this->xEnd_; ++i, ++j) {
284 error = std::fabs(value(*i) - *j);
285 maxError = std::max(maxError, error);
286 }
287 return maxError;
288 }
289
290 private:
291 class XABRError : public CostFunction {
292 public:
293 explicit XABRError(XABRInterpolationImpl *xabr) : xabr_(xabr) {}
294
295 Real value(const Array& x) const override {
296 const Array y = Model().direct(x, xabr_->paramIsFixed_,
298 for (Size i = 0; i < xabr_->params_.size(); ++i)
299 xabr_->params_[i] = y[i];
302 }
303
304 Array values(const Array& x) const override {
305 const Array y = Model().direct(x, xabr_->paramIsFixed_,
307 for (Size i = 0; i < xabr_->params_.size(); ++i)
308 xabr_->params_[i] = y[i];
310 return xabr_->interpolationErrors();
311 }
312
313 private:
315 };
316 ext::shared_ptr<EndCriteria> endCriteria_;
317 ext::shared_ptr<OptimizationMethod> optMethod_;
319 const bool useMaxError_;
324};
325
326} // namespace detail
327} // namespace QuantLib
328
329#endif
Black formula.
1-D array used in linear algebra.
Definition: array.hpp:52
Real * iterator
Definition: array.hpp:123
Size size() const
dimension of the array
Definition: array.hpp:495
Cost function abstract class for optimization problem.
Halton low-discrepancy sequence generator.
Definition: haltonrsg.hpp:43
const sample_type & nextSequence() const
Definition: haltonrsg.cpp:56
basic template implementation
templateImpl(const I1 &xBegin, const I1 &xEnd, const I2 &yBegin, const int requiredPoints=2)
base class for 1-D interpolations.
Levenberg-Marquardt optimization method.
No constraint.
Definition: constraint.hpp:79
template class providing a null value for a given type.
Definition: null.hpp:76
Constrained optimization problem.
Definition: problem.hpp:42
const Array & currentValue() const
current value of the local minimum
Definition: problem.hpp:81
Parameterized cost function.
virtual Array include(const Array &projectedParameters) const
returns whole set of parameters corresponding to the set
Definition: projection.cpp:67
virtual Array project(const Array &parameters) const
returns the subset of free parameters corresponding
Definition: projection.cpp:54
ext::shared_ptr< typename Model::type > modelInstance_
XABRCoeffHolder(const Time t, const Real &forward, const std::vector< Real > &params, const std::vector< bool > &paramIsFixed, std::vector< Real > addParams)
virtual ~XABRCoeffHolder()=default
Real value(const Array &x) const override
method to overload to compute the cost function value in x
Array values(const Array &x) const override
method to overload to compute the cost function values in x
ext::shared_ptr< EndCriteria > endCriteria_
XABRInterpolationImpl(const I1 &xBegin, const I1 &xEnd, const I2 &yBegin, Time t, const Real &forward, const std::vector< Real > &params, const std::vector< bool > &paramIsFixed, bool vegaWeighted, ext::shared_ptr< EndCriteria > endCriteria, ext::shared_ptr< OptimizationMethod > optMethod, const Real errorAccept, const bool useMaxError, const Size maxGuesses, const std::vector< Real > &addParams=std::vector< Real >(), VolatilityType volatilityType=VolatilityType::ShiftedLognormal)
ext::shared_ptr< OptimizationMethod > optMethod_
Abstract constraint class.
output manipulators
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_FAIL(message)
throw an error (possibly with file and line information)
Definition: errors.hpp:92
#define QL_MAX_REAL
Definition: qldefines.hpp:176
#define QL_MIN_REAL
Definition: qldefines.hpp:175
Real Time
continuous quantity with 1-year units
Definition: types.hpp:62
QL_REAL Real
real number
Definition: types.hpp:50
std::size_t Size
size of a container
Definition: types.hpp:58
Halton low-discrepancy sequence generator.
base class for 1-D interpolations
Levenberg-Marquardt optimization method.
Abstract optimization method class.
Definition: any.hpp:35
Matrix inverse(const Matrix &m)
Definition: matrix.cpp:44
STL namespace.
null values
ext::shared_ptr< YieldTermStructure > r
Cost function utility.
Simplex optimization method.
volatility types