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