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