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
20#include <ql/termstructures/volatility/kahalesmilesection.hpp>
21
22using std::sqrt;
23
24namespace QuantLib {
25
26 KahaleSmileSection::KahaleSmileSection(const ext::shared_ptr<SmileSection>& source,
27 const Real atm,
28 const bool interpolate,
29 const bool exponentialExtrapolation,
30 const bool deleteArbitragePoints,
31 const std::vector<Real>& moneynessGrid,
32 const Real gap,
33 const int forcedLeftIndex,
34 const int forcedRightIndex)
35 : SmileSection(*source), source_(source), moneynessGrid_(moneynessGrid), gap_(gap),
36 interpolate_(interpolate), exponentialExtrapolation_(exponentialExtrapolation),
37 forcedLeftIndex_(forcedLeftIndex), forcedRightIndex_(forcedRightIndex) {
38
39 // only shifted lognormal smile sections are supported
40
41 QL_REQUIRE(source->volatilityType() == ShiftedLognormal,
42 "KahaleSmileSection only supports shifted lognormal source sections");
43
44 ssutils_ = ext::make_shared<SmileSectionUtils>(
45 args&: *source, args: moneynessGrid, args: atm, args: deleteArbitragePoints);
46
47 moneynessGrid_ = ssutils_->moneyGrid();
48 k_ = ssutils_->strikeGrid();
49 c_ = ssutils_->callPrices();
50 f_ = ssutils_->atmLevel();
51
52 // for shifted smile sections we shift the forward and the strikes
53 // and do as if we were in a lognormal setting
54
55 for (Real& i : k_) {
56 i += source_->shift();
57 }
58
59 f_ += source_->shift();
60
61 compute();
62 }
63
64 void KahaleSmileSection::compute() {
65
66 std::pair<Size, Size> afIdx = ssutils_->arbitragefreeIndices();
67 leftIndex_ = afIdx.first;
68 rightIndex_ = afIdx.second;
69
70 cFunctions_ = std::vector<ext::shared_ptr<cFunction> >(
71 rightIndex_ - leftIndex_ + 2);
72
73 // extrapolation in the leftmost interval
74
75 Brent brent;
76 bool success;
77 Real secl = 0.0;
78
79 do {
80 success = true;
81 try {
82 Real k1 = k_[leftIndex_];
83 Real c1 = c_[leftIndex_];
84 Real c0 = c_[0];
85 secl = (c_[leftIndex_] - c_[0]) / (k_[leftIndex_] - k_[0]);
86 Real sec = (c_[leftIndex_ + 1] - c_[leftIndex_]) /
87 (k_[leftIndex_ + 1] - k_[leftIndex_]);
88 Real c1p;
89 if (interpolate_)
90 c1p = (secl + sec) / 2;
91 else {
92 c1p = -source_->digitalOptionPrice(strike: k1 - source_->shift() + gap_ / 2.0, type: Option::Call, discount: 1.0, gap: gap_);
93 QL_REQUIRE(secl < c1p && c1p <= 0.0, "dummy");
94 // can not extrapolate so throw exception which is caught
95 // below
96 }
97 sHelper1 sh1(k1, c0, c1, c1p);
98 Real s = brent.solve(f: sh1, QL_KAHALE_ACC, guess: 0.20, xMin: 0.00,
99 QL_KAHALE_SMAX); // numerical parameters
100 // hardcoded here
101 sh1(s);
102 ext::shared_ptr<cFunction> cFct1(
103 new cFunction(sh1.f_, s, 0.0, sh1.b_));
104 cFunctions_[0] = cFct1;
105 // sanity check - in rare cases we can get digitials
106 // which are not monotonic or greater than 1.0
107 // due to numerical effects. Move to the next index in
108 // these cases.
109 Real dig = digitalOptionPrice(strike: (k1 - source_->shift()) / 2.0, type: Option::Call, discount: 1.0, gap: gap_);
110 QL_REQUIRE(dig >= -c1p && dig <= 1.0, "dummy");
111 if(static_cast<int>(leftIndex_) < forcedLeftIndex_) {
112 leftIndex_++;
113 success = false;
114 }
115 }
116 catch (...) {
117 leftIndex_++;
118 success = false;
119 }
120 } while (!success && leftIndex_ < rightIndex_);
121
122 QL_REQUIRE(
123 leftIndex_ < rightIndex_,
124 "can not extrapolate to left, right index of af region reached ("
125 << rightIndex_ << ")");
126
127 // interpolation
128
129 Real cp0 = 0.0, cp1 = 0.0;
130
131 if (interpolate_) {
132
133 for (Size i = leftIndex_; i < rightIndex_; i++) {
134 Real k0 = k_[i];
135 Real k1 = k_[i + 1];
136 Real c0 = c_[i];
137 Real c1 = c_[i + 1];
138 Real sec = (c_[i + 1] - c_[i]) / (k_[i + 1] - k_[i]);
139 if (i == leftIndex_)
140 cp0 = leftIndex_ > 0 ? (secl + sec) / 2.0 : sec;
141 Real secr;
142 if (i == rightIndex_ - 1)
143 secr = 0.0;
144 else
145 secr = (c_[i + 2] - c_[i + 1]) / (k_[i + 2] - k_[i + 1]);
146 cp1 = (sec + secr) / 2.0;
147 aHelper ah(k0, k1, c0, c1, cp0, cp1);
148 Real a;
149 bool valid = false;
150 try {
151 a = brent.solve(
152 f: ah, QL_KAHALE_ACC, guess: 0.5 * (cp1 + (1.0 + cp0)),
153 xMin: cp1 + QL_KAHALE_EPS, xMax: 1.0 + cp0 - QL_KAHALE_EPS);
154 // numerical parameters hardcoded here
155 valid = true;
156 }
157 catch (...) {
158 // delete the right point of the interval where we try to
159 // interpolate
160 moneynessGrid_.erase(position: moneynessGrid_.begin() + (i + 1));
161 k_.erase(position: k_.begin() + (i + 1));
162 c_.erase(position: c_.begin() + (i + 1));
163 cFunctions_.erase(position: cFunctions_.begin() + (i + 1));
164 rightIndex_--;
165 i--;
166 }
167 if (valid) {
168 ah(a);
169 ext::shared_ptr<cFunction> cFct(
170 new cFunction(ah.f_, ah.s_, a, ah.b_));
171 cFunctions_[leftIndex_ > 0 ? i - leftIndex_ + 1 : 0] = cFct;
172 cp0 = cp1;
173 }
174 }
175 }
176
177 // extrapolation of right wing
178
179 do {
180 success = true;
181 try {
182 Real k0 = k_[rightIndex_];
183 Real c0 = c_[rightIndex_];
184 Real cp0;
185 if (interpolate_)
186 cp0 = 0.5 * (c_[rightIndex_] - c_[rightIndex_ - 1]) /
187 (k_[rightIndex_] - k_[rightIndex_ - 1]);
188 else {
189 cp0 = -source_->digitalOptionPrice(
190 strike: k0 - shift() - gap_ / 2.0, type: Option::Call, discount: 1.0, gap: gap_);
191 }
192 ext::shared_ptr<cFunction> cFct;
193 if (exponentialExtrapolation_) {
194 QL_REQUIRE(-cp0 / c0 > 0.0, "dummy"); // this is caught
195 // below
196 cFct = ext::make_shared<cFunction>(
197 args: -cp0 / c0, args: std::log(x: c0) - cp0 / c0 * k0);
198 } else {
199 sHelper sh(k0, c0, cp0);
200 Real s;
201 s = brent.solve(f: sh, QL_KAHALE_ACC, guess: 0.20, xMin: 0.0,
202 QL_KAHALE_SMAX); // numerical parameters
203 // hardcoded here
204 sh(s);
205 cFct = ext::make_shared<cFunction>(
206 args&: sh.f_, args&: s, args: 0.0, args: 0.0);
207 }
208 cFunctions_[rightIndex_ - leftIndex_ + 1] = cFct;
209 }
210 catch (...) {
211 rightIndex_--;
212 success = false;
213 }
214 if(static_cast<int>(rightIndex_) > forcedRightIndex_) {
215 rightIndex_--;
216 success = false;
217 }
218 } while (!success && rightIndex_ > leftIndex_);
219
220 QL_REQUIRE(
221 leftIndex_ < rightIndex_,
222 "can not extrapolate to right, left index of af region reached ("
223 << leftIndex_ << ")");
224 }
225
226 Real KahaleSmileSection::optionPrice(Rate strike, Option::Type type,
227 Real discount) const {
228 // option prices are directly available, so implement this function
229 // rather than use smileSection
230 // standard implementation
231 Real shifted_strike = std::max(a: strike + shift(), QL_KAHALE_EPS);
232 int i = index(strike: shifted_strike);
233 if (interpolate_ ||
234 (i == 0 || i == (int)(rightIndex_ - leftIndex_ + 1)))
235 return discount *
236 (type == Option::Call
237 ? (*cFunctions_[i])(shifted_strike)
238 : (*cFunctions_[i])(shifted_strike) + shifted_strike - f_);
239 else
240 return source_->optionPrice(strike, type, discount);
241 }
242
243 Real KahaleSmileSection::volatilityImpl(Rate strike) const {
244 Real shifted_strike = std::max(a: strike + shift(), QL_KAHALE_EPS);
245 int i = index(strike: shifted_strike);
246 if (!interpolate_ &&
247 !(i == 0 || i == (int)(rightIndex_ - leftIndex_ + 1)))
248 return source_->volatility(strike);
249 Real c = (*cFunctions_[i])(shifted_strike);
250 Real vol = 0.0;
251 try {
252 Option::Type type = shifted_strike >= f_ ? Option::Call : Option::Put;
253 vol = blackFormulaImpliedStdDev(
254 optionType: type, strike: shifted_strike, forward: f_,
255 blackPrice: type == Option::Put ? strike - f_ + c : c) /
256 sqrt(x: exerciseTime());
257 }
258 catch (...) {
259 }
260 return vol;
261 }
262
263 Size KahaleSmileSection::index(Rate strike) const {
264 int i =
265 static_cast<int>(std::upper_bound(first: k_.begin(), last: k_.end(), val: strike) -
266 k_.begin()) -
267 static_cast<int>(leftIndex_);
268 return std::max(
269 a: std::min(a: i, b: static_cast<int>(rightIndex_ - leftIndex_ + 1)), b: 0);
270 }
271}
272

source code of quantlib/ql/termstructures/volatility/kahalesmilesection.cpp