1 /* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ 2 3 /* 4 Copyright (C) 2014 Master IMAFA - Polytech'Nice Sophia - Université de Nice Sophia Antipolis 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/experimental/exoticoptions/analyticholderextensibleoptionengine.hpp> 21 #include <ql/math/distributions/bivariatenormaldistribution.hpp> 22 #include <ql/exercise.hpp> 23 24 using std::pow; 25 using std::log; 26 using std::exp; 27 using std::sqrt; 28 29 namespace QuantLib { 30 AnalyticHolderExtensibleOptionEngine(const ext::shared_ptr<GeneralizedBlackScholesProcess> & process)31 AnalyticHolderExtensibleOptionEngine::AnalyticHolderExtensibleOptionEngine( 32 const ext::shared_ptr<GeneralizedBlackScholesProcess>& process) 33 : process_(process) { 34 registerWith(process_); 35 } 36 calculate() const37 void AnalyticHolderExtensibleOptionEngine::calculate() const { 38 //Spot 39 Real S = process_->x0(); 40 Real r = riskFreeRate(); 41 Real b = r - dividendYield(); 42 Real X1 = strike(); 43 Real X2 = arguments_.secondStrike; 44 Time T2 = secondExpiryTime(); 45 Time t1 = firstExpiryTime(); 46 Real A = arguments_.premium; 47 48 49 Real z1 = this->z1(); 50 51 Real z2 = this->z2(); 52 53 Real rho = sqrt(t1 / T2); 54 55 56 ext::shared_ptr<PlainVanillaPayoff> payoff = 57 ext::dynamic_pointer_cast<PlainVanillaPayoff>(arguments_.payoff); 58 59 //QuantLib requires sigma * sqrt(T) rather than just sigma/volatility 60 Real vol = volatility(); 61 62 //calculate dividend discount factor assuming continuous compounding (e^-rt) 63 DiscountFactor growth = dividendDiscount(t1); 64 //calculate payoff discount factor assuming continuous compounding 65 DiscountFactor discount = riskFreeDiscount(t1); 66 Real result = 0; 67 Real minusInf=-std::numeric_limits<Real>::infinity(); 68 69 Real y1 = this->y1(payoff->optionType()), 70 y2 = this->y2(payoff->optionType()); 71 if (payoff->optionType() == Option::Call) { 72 //instantiate payoff function for a call 73 ext::shared_ptr<PlainVanillaPayoff> vanillaCallPayoff = 74 ext::make_shared<PlainVanillaPayoff>(Option::Call, X1); 75 Real BSM = BlackScholesCalculator(vanillaCallPayoff, S, growth, vol*sqrt(t1), discount).value(); 76 result = BSM 77 + S*exp((b - r)*T2)*M2(y1, y2, minusInf, z1, rho) 78 - X2*exp(-r*T2)*M2(y1 - vol*sqrt(t1), y2 - vol*sqrt(t1), minusInf, z1 - vol*sqrt(T2), rho) 79 - S*exp((b - r)*t1)*N2(y1, z2) + X1*exp(-r*t1)*N2(y1 - vol*sqrt(t1), z2 - vol*sqrt(t1)) 80 - A*exp(-r*t1)*N2(y1 - vol*sqrt(t1), y2 - vol*sqrt(t1)); 81 } else { 82 //instantiate payoff function for a call 83 ext::shared_ptr<PlainVanillaPayoff> vanillaPutPayoff = 84 ext::make_shared<PlainVanillaPayoff>(Option::Put, X1); 85 result = BlackScholesCalculator(vanillaPutPayoff, S, growth, vol*sqrt(t1), discount).value() 86 - S*exp((b - r)*T2)*M2(y1, y2, minusInf, -z1, rho) 87 + X2*exp(-r*T2)*M2(y1 - vol*sqrt(t1), y2 - vol*sqrt(t1), minusInf, -z1 + vol*sqrt(T2), rho) 88 + S*exp((b - r)*t1)*N2(z2, y2) - X1*exp(-r*t1)*N2(z2 - vol*sqrt(t1), y2 - vol*sqrt(t1)) 89 - A*exp(-r*t1)*N2(y1 - vol*sqrt(t1), y2 - vol*sqrt(t1)); 90 } 91 this->results_.value = result; 92 } 93 I1Call() const94 Real AnalyticHolderExtensibleOptionEngine::I1Call() const { 95 Real Sv = process_->x0(); 96 Real A = arguments_.premium; 97 98 if(A==0) 99 { 100 return 0; 101 } 102 else 103 { 104 BlackScholesCalculator bs = bsCalculator(Sv, Option::Call); 105 Real ci = bs.value(); 106 Real dc = bs.delta(); 107 108 Real yi = ci - A; 109 //da/ds = 0 110 Real di = dc - 0; 111 Real epsilon = 0.001; 112 113 //Newton-Raphson process 114 while (std::fabs(yi) > epsilon){ 115 Sv = Sv - yi / di; 116 117 bs = bsCalculator(Sv, Option::Call); 118 ci = bs.value(); 119 dc = bs.delta(); 120 121 yi = ci - A; 122 di = dc - 0; 123 } 124 return Sv; 125 } 126 } 127 I2Call() const128 Real AnalyticHolderExtensibleOptionEngine::I2Call() const { 129 Real Sv = process_->x0(); 130 Real X1 = strike(); 131 Real X2 = arguments_.secondStrike; 132 Real A = arguments_.premium; 133 Time T2 = secondExpiryTime(); 134 Time t1 = firstExpiryTime(); 135 Real r=riskFreeRate(); 136 137 Real val=X1-X2*std::exp(-r*(T2-t1)); 138 if(A< val){ 139 return std::numeric_limits<Real>::infinity(); 140 } else { 141 BlackScholesCalculator bs = bsCalculator(Sv, Option::Call); 142 Real ci = bs.value(); 143 Real dc = bs.delta(); 144 145 Real yi = ci - A - Sv + X1; 146 //da/ds = 1 147 Real di = dc - 1; 148 Real epsilon = 0.001; 149 150 //Newton-Raphson process 151 while (std::fabs(yi) > epsilon){ 152 Sv = Sv - yi / di; 153 154 bs = bsCalculator(Sv, Option::Call); 155 ci = bs.value(); 156 dc = bs.delta(); 157 158 yi = ci - A - Sv + X1; 159 di = dc - 1; 160 } 161 return Sv; 162 } 163 } 164 I1Put() const165 Real AnalyticHolderExtensibleOptionEngine::I1Put() const { 166 Real Sv = process_->x0(); 167 //Srtike 168 Real X1 = strike(); 169 //Premium 170 Real A = arguments_.premium; 171 172 BlackScholesCalculator bs = bsCalculator(Sv, Option::Put); 173 Real pi = bs.value(); 174 Real dc = bs.delta(); 175 176 Real yi = pi - A + Sv - X1; 177 //da/ds = 1 178 Real di = dc - 1; 179 Real epsilon = 0.001; 180 181 //Newton-Raphson prosess 182 while (std::fabs(yi) > epsilon){ 183 Sv = Sv - yi / di; 184 185 bs = bsCalculator(Sv, Option::Put); 186 pi = bs.value(); 187 dc = bs.delta(); 188 189 yi = pi - A + Sv - X1; 190 di = dc - 1; 191 } 192 return Sv; 193 } 194 I2Put() const195 Real AnalyticHolderExtensibleOptionEngine::I2Put() const { 196 Real Sv = process_->x0(); 197 Real A = arguments_.premium; 198 if(A==0){ 199 return std::numeric_limits<Real>::infinity(); 200 } 201 else{ 202 BlackScholesCalculator bs = bsCalculator(Sv, Option::Put); 203 Real pi = bs.value(); 204 Real dc = bs.delta(); 205 206 Real yi = pi - A; 207 //da/ds = 0 208 Real di = dc - 0; 209 Real epsilon = 0.001; 210 211 //Newton-Raphson prosess 212 while (std::fabs(yi) > epsilon){ 213 Sv = Sv - yi / di; 214 215 bs = bsCalculator(Sv, Option::Put); 216 pi = bs.value(); 217 dc = bs.delta(); 218 219 yi = pi - A; 220 di = dc - 0; 221 } 222 return Sv; 223 } 224 } 225 226 bsCalculator(Real spot,Option::Type optionType) const227 BlackScholesCalculator AnalyticHolderExtensibleOptionEngine::bsCalculator( 228 Real spot, Option::Type optionType) const { 229 //Real spot = process_->x0(); 230 Real vol; 231 DiscountFactor growth; 232 DiscountFactor discount; 233 Real X2 = arguments_.secondStrike; 234 Time T2 = secondExpiryTime(); 235 Time t1 = firstExpiryTime(); 236 Time t = T2 - t1; 237 238 //payoff 239 ext::shared_ptr<PlainVanillaPayoff > vanillaPayoff = 240 ext::make_shared<PlainVanillaPayoff>(optionType, X2); 241 242 //QuantLib requires sigma * sqrt(T) rather than just sigma/volatility 243 vol = volatility() * std::sqrt(t); 244 //calculate dividend discount factor assuming continuous compounding (e^-rt) 245 growth = dividendDiscount(t); 246 //calculate payoff discount factor assuming continuous compounding 247 discount = riskFreeDiscount(t); 248 249 BlackScholesCalculator bs(vanillaPayoff, spot, growth, vol, discount); 250 return bs; 251 } 252 M2(Real a,Real b,Real c,Real d,Real rho) const253 Real AnalyticHolderExtensibleOptionEngine::M2(Real a, Real b, Real c, Real d, Real rho) const { 254 BivariateCumulativeNormalDistributionDr78 CmlNormDist(rho); 255 return CmlNormDist(b, d) - CmlNormDist(a, d) - CmlNormDist(b, c) + CmlNormDist(a,c); 256 } 257 N2(Real a,Real b) const258 Real AnalyticHolderExtensibleOptionEngine::N2(Real a, Real b) const { 259 CumulativeNormalDistribution NormDist; 260 return NormDist(b) - NormDist(a); 261 } 262 strike() const263 Real AnalyticHolderExtensibleOptionEngine::strike() const { 264 ext::shared_ptr<PlainVanillaPayoff> payoff = 265 ext::dynamic_pointer_cast<PlainVanillaPayoff>(arguments_.payoff); 266 QL_REQUIRE(payoff, "non-plain payoff given"); 267 return payoff->strike(); 268 } 269 firstExpiryTime() const270 Time AnalyticHolderExtensibleOptionEngine::firstExpiryTime() const { 271 return process_->time(arguments_.exercise->lastDate()); 272 } 273 secondExpiryTime() const274 Time AnalyticHolderExtensibleOptionEngine::secondExpiryTime() const { 275 return process_->time(arguments_.secondExpiryDate); 276 } 277 volatility() const278 Volatility AnalyticHolderExtensibleOptionEngine::volatility() const { 279 return process_->blackVolatility()->blackVol(firstExpiryTime(), strike()); 280 } riskFreeRate() const281 Rate AnalyticHolderExtensibleOptionEngine::riskFreeRate() const { 282 return process_->riskFreeRate()->zeroRate(firstExpiryTime(), Continuous, 283 NoFrequency); 284 } dividendYield() const285 Rate AnalyticHolderExtensibleOptionEngine::dividendYield() const { 286 return process_->dividendYield()->zeroRate(firstExpiryTime(), 287 Continuous, NoFrequency); 288 } 289 dividendDiscount(Time t) const290 DiscountFactor AnalyticHolderExtensibleOptionEngine::dividendDiscount(Time t) const { 291 return process_->dividendYield()->discount(t); 292 } 293 riskFreeDiscount(Time t) const294 DiscountFactor AnalyticHolderExtensibleOptionEngine::riskFreeDiscount(Time t) const { 295 return process_->riskFreeRate()->discount(t); 296 } 297 y1(Option::Type type) const298 Real AnalyticHolderExtensibleOptionEngine::y1(Option::Type type) const { 299 Real S = process_->x0(); 300 Real I2 = (type == Option::Call) ? I2Call() : I2Put(); 301 302 Real b = riskFreeRate() - dividendYield(); 303 Real vol = volatility(); 304 Time t1 = firstExpiryTime(); 305 306 return (log(S / I2) + (b + pow(vol, 2) / 2)*t1) / (vol*sqrt(t1)); 307 } 308 y2(Option::Type type) const309 Real AnalyticHolderExtensibleOptionEngine::y2(Option::Type type) const { 310 Real S = process_->x0(); 311 Real I1 = (type == Option::Call) ? I1Call() : I1Put(); 312 313 Real b = riskFreeRate() - dividendYield(); 314 Real vol = volatility(); 315 Time t1 = firstExpiryTime(); 316 317 return (log(S / I1) + (b + pow(vol, 2) / 2)*t1) / (vol*sqrt(t1)); 318 } 319 z1() const320 Real AnalyticHolderExtensibleOptionEngine::z1() const { 321 Real S = process_->x0(); 322 Real X2 = arguments_.secondStrike; 323 Real b = riskFreeRate() - dividendYield(); 324 Real vol = volatility(); 325 Time T2 = secondExpiryTime(); 326 327 return (log(S / X2) + (b + pow(vol, 2) / 2)*T2) / (vol*sqrt(T2)); 328 } 329 z2() const330 Real AnalyticHolderExtensibleOptionEngine::z2() const { 331 Real S = process_->x0(); 332 Real X1 = strike(); 333 334 Real b = riskFreeRate() - dividendYield(); 335 Real vol = volatility(); 336 Time t1 = firstExpiryTime(); 337 338 return (log(S / X1) + (b + pow(vol, 2) / 2)*t1) / (vol*sqrt(t1)); 339 } 340 341 } 342