45 struct close_enough_to {
48 explicit close_enough_to(
Real y,
Size n=42) :
y(
y),
n(
n) {}
54 class AndreasenHugeCostFunction :
public CostFunction {
56 AndreasenHugeCostFunction(
59 Array lnMarketStrikes,
61 const ext::shared_ptr<FdmMesherComposite>& mesher,
64 : marketNPVs_(
std::move(marketNPVs)), marketVegas_(
std::move(marketVegas)),
65 lnMarketStrikes_(
std::move(lnMarketStrikes)), previousNPVs_(
std::move(previousNPVs)),
66 mesher_(mesher), nGridPoints_(mesher->layout()->size()), dT_(dT),
67 interpolationType_((lnMarketStrikes_.size() > 1) ?
69 AndreasenHugeVolatilityInterpl::PiecewiseConstant),
70 dxMap_(FirstDerivativeOp(0,
mesher_)), dxxMap_(SecondDerivativeOp(0,
mesher_)),
71 d2CdK2_(dxMap_.mult(Array(mesher->layout()->size(), -1.0)).add(dxxMap_)),
74 Array d2CdK2(
const Array& c)
const {
75 return d2CdK2_.apply(c);
78 Array solveFor(
Time dT,
const Array& sig,
const Array&
b)
const {
80 Array x(lnMarketStrikes_.size());
81 Interpolation sigInterpl;
83 switch (interpolationType_) {
85 sigInterpl = CubicNaturalSpline(
86 lnMarketStrikes_.begin(), lnMarketStrikes_.end(),
90 sigInterpl = LinearInterpolation(
91 lnMarketStrikes_.begin(), lnMarketStrikes_.end(),
95 for (
Size i=0; i < x.size()-1; ++i)
96 x[i] = 0.5*(lnMarketStrikes_[i] + lnMarketStrikes_[i+1]);
97 x.back() = lnMarketStrikes_.back();
99 sigInterpl = BackwardFlatInterpolation(
100 x.begin(), x.end(), sig.begin());
103 QL_FAIL(
"unknown interpolation type");
106 Array z(
mesher_->layout()->size());
107 for (
const auto& iter : *
mesher_->layout()) {
108 const Size i = iter.index();
111 const Real vol = sigInterpl(
112 std::min(std::max(lnStrike, lnMarketStrikes_.front()),
113 lnMarketStrikes_.back()),
true);
118 mapT_.axpyb(z, dxMap_, dxxMap_.mult(-z), Array());
119 return mapT_.mult(Array(z.size(), dT)).solve_splitting(
b, 1.0);
123 Array
apply(
const Array& c)
const {
124 return -mapT_.apply(c);
127 Array values(
const Array& sig)
const override {
128 Array newNPVs = solveFor(dT_, sig, previousNPVs_);
130 const std::vector<Real>& gridPoints =
131 mesher_->getFdm1dMeshers().front()->locations();
133 const MonotonicCubicNaturalSpline interpl(
134 gridPoints.begin(), gridPoints.end(), newNPVs.begin());
136 Array retVal(lnMarketStrikes_.size());
137 for (
Size i=0; i < retVal.size(); ++i) {
138 const Real strike = lnMarketStrikes_[i];
139 retVal[i] = interpl(strike) - marketNPVs_[i];
144 Array vegaCalibrationError(
const Array& sig)
const {
145 return values(sig)/marketVegas_;
148 Array initialValues()
const {
149 return Array(lnMarketStrikes_.size(), 0.25);
154 const Array marketNPVs_, marketVegas_;
155 const Array lnMarketStrikes_, previousNPVs_;
156 const ext::shared_ptr<FdmMesherComposite>
mesher_;
157 const Size nGridPoints_;
162 const FirstDerivativeOp dxMap_;
163 const TripleBandLinearOp dxxMap_;
164 const TripleBandLinearOp d2CdK2_;
165 mutable TripleBandLinearOp mapT_;
168 class CombinedCostFunction :
public CostFunction {
170 CombinedCostFunction(ext::shared_ptr<AndreasenHugeCostFunction> putCostFct,
171 ext::shared_ptr<AndreasenHugeCostFunction> callCostFct)
172 : putCostFct_(
std::move(putCostFct)), callCostFct_(
std::move(callCostFct)) {}
174 Array values(
const Array& sig)
const override {
175 if ((putCostFct_ !=
nullptr) && (callCostFct_ !=
nullptr)) {
176 const Array pv = putCostFct_->values(sig);
177 const Array cv = callCostFct_->values(sig);
179 Array retVal(pv.size() + cv.size());
180 std::copy(pv.begin(), pv.end(), retVal.begin());
181 std::copy(cv.begin(), cv.end(), retVal.begin() + cv.size());
184 }
else if (putCostFct_ !=
nullptr)
185 return putCostFct_->values(sig);
186 else if (callCostFct_ !=
nullptr)
187 return callCostFct_->values(sig);
189 QL_FAIL(
"internal error: cost function not set");
192 Array initialValues()
const {
193 if ((putCostFct_ !=
nullptr) && (callCostFct_ !=
nullptr))
194 return 0.5*( putCostFct_->initialValues()
195 + callCostFct_->initialValues());
196 else if (putCostFct_ !=
nullptr)
197 return putCostFct_->initialValues();
198 else if (callCostFct_ !=
nullptr)
199 return callCostFct_->initialValues();
201 QL_FAIL(
"internal error: cost function not set");
205 const ext::shared_ptr<AndreasenHugeCostFunction> putCostFct_;
206 const ext::shared_ptr<AndreasenHugeCostFunction> callCostFct_;
220 ext::shared_ptr<OptimizationMethod> optimizationMethod,
222 : spot_(
std::move(spot)),
rTS_(
std::move(rTS)), qTS_(
std::move(qTS)),
223 interpolationType_(interplationType), calibrationType_(calibrationType),
224 nGridPoints_(nGridPoints), minStrike_(_minStrike), maxStrike_(_maxStrike),
225 optimizationMethod_(
std::move(optimizationMethod)), endCriteria_(endCriteria) {
226 QL_REQUIRE(nGridPoints > 2 && !calibrationSet.empty(),
"undefined grid or calibration set");
228 std::set<Real> strikes;
229 std::set<Date> expiries;
232 for (
const auto& i : calibrationSet) {
234 const ext::shared_ptr<Exercise> exercise = i.first->exercise();
237 "European option required");
239 const Date expiry = exercise->lastDate();
240 expiries.insert(expiry);
242 const ext::shared_ptr<PlainVanillaPayoff>
payoff =
243 ext::dynamic_pointer_cast<PlainVanillaPayoff>(i.first->payoff());
248 strikes.insert(strike);
255 strikes_.assign(strikes.begin(), strikes.end());
256 expiries_.assign(expiries.begin(), expiries.end());
262 expiries.size(), std::vector<Size>(strikes.size(),
Null<Size>()));
264 for (
Size i=0; i < calibrationSet.size(); ++i) {
266 calibrationSet[i].first->exercise()->lastDate();
268 const Size l = std::distance(expiries.begin(), expiries.lower_bound(expiry));
271 ext::dynamic_pointer_cast<PlainVanillaPayoff>(
272 calibrationSet[i].first->payoff())->strike();
276 close_enough_to(strike)));
286 ext::shared_ptr<AndreasenHugeCostFunction>
289 const Array& previousNPVs)
const {
294 return ext::shared_ptr<AndreasenHugeCostFunction>();
302 const Size nOptions = std::count_if(
305 [=](
Size n){
return n != null; });
307 Array lnMarketStrikes(nOptions),
308 marketNPVs(nOptions), marketVegas(nOptions);
317 const Real stdDev = vol*std::sqrt(expiryTime);
323 const Real vega = calculator.
vega(expiryTime);
325 marketNPVs[k] = npv/(discount*
fwd);
326 marketVegas[k] = vega/(discount*
fwd);
331 return ext::make_shared<AndreasenHugeCostFunction>(
344 "max strike must be greater than min strike");
354 ext::make_shared<FdmMesherComposite>(
355 ext::make_shared<Concentrating1dMesher>(
359 std::pair<Real, Real>(0.0, 0.025)));
368 minError_ = std::numeric_limits<Real>::max();
383 const ext::shared_ptr<AndreasenHugeCostFunction> putCostFct =
385 const ext::shared_ptr<AndreasenHugeCostFunction> callCostFct =
388 CombinedCostFunction costFunction(putCostFct, callCostFct);
392 positiveConstraint, costFunction.initialValues());
399 npvPuts, npvCalls, sig,
408 const Array vegaPutDiffs =
409 putCostFct->vegaCalibrationError(sig);
410 const Array vegaCallDiffs =
411 callCostFct->vegaCalibrationError(sig);
416 for (
Size j=0; j < vegaDiffs.
size(); ++j)
417 vegaDiffs[j] = std::fabs(
422 vegaDiffs =
Abs(putCostFct->vegaCalibrationError(sig));
425 vegaDiffs =
Abs(callCostFct->vegaCalibrationError(sig));
428 QL_FAIL(
"unknown calibration type");
432 std::accumulate(vegaDiffs.
begin(), vegaDiffs.
end(),
Real(0.0));
434 *std::min_element(vegaDiffs.
begin(), vegaDiffs.
end()));
436 *std::max_element(vegaDiffs.
begin(), vegaDiffs.
end()));
438 if (putCostFct !=
nullptr)
439 npvPuts = putCostFct->solveFor(
dT_[i], sig, npvPuts);
440 if (callCostFct !=
nullptr)
441 npvCalls= callCostFct->solveFor(
dT_[i], sig, npvCalls);
470 ext::tuple<Real, Real, Real>
485 Real strike,
const TimeValueCacheType::const_iterator&
f)
const {
487 const Real fwd = ext::get<0>(
f->second);
488 const Real k = std::log(strike /
fwd);
493 return (*(ext::get<2>(
f->second)))(
s);
516 const Real fwd = ext::get<0>(
f->second);
522 price = price + strike/
fwd - 1.0;
524 price = 1.0 - strike/
fwd + price;
532 ext::shared_ptr<Array> prices(
544 QL_FAIL(
"unknown calibration type");
551 ext::make_shared<CubicNaturalSpline>(
563 const Array& previousNPVs =
567 const ext::shared_ptr<AndreasenHugeCostFunction> costFunction
573 const Array cAtJ = costFunction->solveFor(dt, sig, previousNPVs);
576 costFunction->solveFor(dt, sig,
578 costFunction->solveFor(dt, sig, previousNPVs)));
580 const Array d2CdK2 = costFunction->d2CdK2(cAtJ);
610 (
gridPoints_[i] > 0.0)? callLocalVol[i] : putLocalVol[i];
620 QL_FAIL(
"unknown calibration type");
627 ext::make_shared<LinearInterpolation>(
Andreasen-Huge local volatility calibration and interpolation.
1-D array used in linear algebra.
backward-flat interpolation between discrete points
Black-formula calculator class.
const Handle< Quote > spot_
ext::tuple< Real, Real, Real > calibrationError() const
void performCalculations() const override
const Handle< YieldTermStructure > qTS_
Array getPriceSlice(Time t, Option::Type optionType) const
std::vector< Real > strikes_
Volatility localVol(Time t, Real strike) const
ext::shared_ptr< FdmMesherComposite > mesher_
Array getLocalVolSlice(Time t, Option::Type optionType) const
TimeValueCacheType priceCache_
const EndCriteria endCriteria_
std::vector< std::pair< ext::shared_ptr< VanillaOption >, ext::shared_ptr< Quote > > > CalibrationSet
std::vector< Date > expiries_
Real getCacheValue(Real strike, const TimeValueCacheType::const_iterator &f) const
AndreasenHugeVolatilityInterpl(const CalibrationSet &calibrationSet, Handle< Quote > spot, Handle< YieldTermStructure > rTS, Handle< YieldTermStructure > qTS, InterpolationType interpolationType=CubicSpline, CalibrationType calibrationType=Call, Size nGridPoints=500, Real minStrike=Null< Real >(), Real maxStrike=Null< Real >(), ext::shared_ptr< OptimizationMethod > optimizationMethod=ext::shared_ptr< OptimizationMethod >(new LevenbergMarquardt), const EndCriteria &endCriteria=EndCriteria(500, 100, 1e-12, 1e-10, 1e-10))
TimeValueCacheType localVolCache_
ext::shared_ptr< AndreasenHugeCostFunction > buildCostFunction(Size iExpiry, Option::Type optionType, const Array &previousNPVs) const
const InterpolationType interpolationType_
std::vector< SingleStepCalibrationResult > calibrationResults_
const CalibrationType calibrationType_
std::vector< Time > expiryTimes_
const ext::shared_ptr< OptimizationMethod > optimizationMethod_
std::vector< std::vector< Size > > calibrationMatrix_
CalibrationSet calibrationSet_
const Handle< YieldTermStructure > rTS_
Size getExerciseTimeIdx(Time t) const
Real optionPrice(Time t, Real strike, Option::Type optionType) const
const Handle< YieldTermStructure > & riskFreeRate() const
1-D array used in linear algebra.
const_iterator end() const
Size size() const
dimension of the array
const_iterator begin() const
Black 1976 calculator class.
Real vega(Time maturity) const
Time yearFraction(const Date &, const Date &, const Date &refPeriodStart=Date(), const Date &refPeriodEnd=Date()) const
Returns the period between two dates as a fraction of year.
Criteria to end optimization process:
Shared handle to an observable.
virtual void calculate() const
template class providing a null value for a given type.
std::pair< iterator, bool > registerWith(const ext::shared_ptr< Observable > &)
Constraint imposing positivity to all arguments
Constrained optimization problem.
const Array & currentValue() const
current value of the local minimum
floating-point comparisons
One-dimensional grid mesher concentrating around critical points.
cubic interpolation between discrete points
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
#define QL_FAIL(message)
throw an error (possibly with file and line information)
Option exercise classes and payoff function.
ext::function< Real(Real)> b
const ext::shared_ptr< YieldTermStructure > rTS_
memory layout of a fdm linear operator
FdmMesher which is a composite of Fdm1dMesher.
const ext::shared_ptr< FdmMesher > mesher_
first derivative linear operator
Real Time
continuous quantity with 1-year units
Real DiscountFactor
discount factor between dates
Real Volatility
volatility
std::size_t Size
size of a container
ext::shared_ptr< QuantLib::Payoff > payoff
Array apply(const Array &x, const ext::function< Real(Real)> &f)
Array Abs(const Array &v)
Array Exp(const Array &v)
Array Sqrt(const Array &v)
bool close_enough(const Quantity &m1, const Quantity &m2, Size n)
second derivative operator
Vanilla option on a single asset.
Interest-rate term structure.