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/analyticpartialtimebarrieroptionengine.hpp> 21 #include <ql/pricingengines/vanilla/analyticeuropeanengine.hpp> 22 #include <ql/math/distributions/bivariatenormaldistribution.hpp> 23 #include <ql/exercise.hpp> 24 25 namespace QuantLib { 26 AnalyticPartialTimeBarrierOptionEngine(const ext::shared_ptr<GeneralizedBlackScholesProcess> & process)27 AnalyticPartialTimeBarrierOptionEngine::AnalyticPartialTimeBarrierOptionEngine( 28 const ext::shared_ptr<GeneralizedBlackScholesProcess>& process) 29 : process_(process) { 30 registerWith(process_); 31 } 32 calculate() const33 void AnalyticPartialTimeBarrierOptionEngine::calculate() const { 34 ext::shared_ptr<PlainVanillaPayoff> payoff = 35 ext::dynamic_pointer_cast<PlainVanillaPayoff>(arguments_.payoff); 36 QL_REQUIRE(payoff, "non-plain payoff given"); 37 QL_REQUIRE(payoff->strike()>0.0, 38 "strike must be positive"); 39 40 Real spot = process_->x0(); 41 QL_REQUIRE(spot >= 0.0, "negative or null underlying given"); 42 43 PartialBarrier::Type barrierType = arguments_.barrierType; 44 PartialBarrier::Range barrierRange = arguments_.barrierRange; 45 46 switch (payoff->optionType()) { 47 //Call Option 48 case Option::Call: 49 switch (barrierType) { 50 case PartialBarrier::DownOut: 51 switch (barrierRange) { 52 case PartialBarrier::Start: 53 results_.value = CA(1); 54 break; 55 case PartialBarrier::EndB1: 56 results_.value = CoB1(); 57 break; 58 case PartialBarrier::EndB2: 59 results_.value = CoB2(PartialBarrier::DownOut); 60 break; 61 default: 62 QL_FAIL("invalid barrier range"); 63 } 64 break; 65 66 case PartialBarrier::DownIn: 67 switch (barrierRange) { 68 case PartialBarrier::Start: 69 results_.value = CIA(1); 70 break; 71 case PartialBarrier::End: 72 QL_FAIL("Down-and-in partial-time end barrier is not implemented"); 73 default: 74 QL_FAIL("invalid barrier range"); 75 } 76 break; 77 78 case PartialBarrier::UpOut: 79 switch (barrierRange) { 80 case PartialBarrier::Start: 81 results_.value = CA(-1); 82 break; 83 case PartialBarrier::EndB1: 84 results_.value = CoB1(); 85 break; 86 case PartialBarrier::EndB2: 87 results_.value = CoB2(PartialBarrier::UpOut); 88 break; 89 default: 90 QL_FAIL("invalid barrier range"); 91 } 92 break; 93 94 case PartialBarrier::UpIn: 95 switch (barrierRange) { 96 case PartialBarrier::Start: 97 results_.value = CIA(-1); 98 break; 99 case PartialBarrier::End: 100 QL_FAIL("Up-and-in partial-time end barrier is not implemented"); 101 default: 102 QL_FAIL("invalid barrier range"); 103 } 104 break; 105 default: 106 QL_FAIL("unknown barrier type"); 107 } 108 break; 109 110 case Option::Put: 111 QL_FAIL("Partial-time barrier Put option is not implemented"); 112 113 default: 114 QL_FAIL("unknown option type"); 115 } 116 } 117 CoB2(PartialBarrier::Type barrierType) const118 Real AnalyticPartialTimeBarrierOptionEngine::CoB2( 119 PartialBarrier::Type barrierType) const { 120 Real result = 0.0; 121 Real b = riskFreeRate()-dividendYield(); 122 if (strike()<barrier()){ 123 switch (barrierType) { 124 case PartialBarrier::DownOut: 125 result = underlying()*std::exp((b-riskFreeRate())*residualTime()); 126 result *= (M(g1(),e1(),rho())-HS(underlying(),barrier(),2*(mu()+1))*M(g3(),-e3(),-rho())); 127 result -= strike()*std::exp(-riskFreeRate()*residualTime())*(M(g2(),e2(),rho())-HS(underlying(),barrier(),2*mu())*M(g4(),-e4(),-rho())); 128 return result; 129 130 case PartialBarrier::UpOut: 131 result = underlying()*std::exp((b-riskFreeRate())*residualTime()); 132 result *= (M(-g1(),-e1(),rho())-HS(underlying(),barrier(),2*(mu()+1))*M(-g3(),e3(),-rho())); 133 result -= strike()*std::exp(-riskFreeRate()*residualTime())*(M(-g2(),-e2(),rho())-HS(underlying(),barrier(),2*mu())*M(-g4(),e4(),-rho())); 134 result -= underlying()*std::exp((b-riskFreeRate())*residualTime())*(M(-d1(),-e1(),rho())-HS(underlying(),barrier(),2*(mu()+1))*M(e3(),-f1(),-rho())); 135 result += strike()*std::exp(-riskFreeRate()*residualTime())*(M(-d2(),-e2(),rho())-HS(underlying(),barrier(),2*mu())*M(e4(),-f2(),-rho())); 136 return result; 137 138 default: 139 QL_FAIL("invalid barrier type"); 140 } 141 } else { 142 QL_FAIL("case of strike>barrier is not implemented for OutEnd B2 type"); 143 } 144 } 145 CoB1() const146 Real AnalyticPartialTimeBarrierOptionEngine::CoB1() const { 147 Real result = 0.0; 148 Real b = riskFreeRate()-dividendYield(); 149 if (strike()>barrier()) { 150 result = underlying()*std::exp((b-riskFreeRate())*residualTime()); 151 result *= (M(d1(),e1(),rho())-HS(underlying(),barrier(),2*(mu()+1))*M(f1(),-e3(),-rho())); 152 result -= (strike()*std::exp(-riskFreeRate()*residualTime()))*(M(d2(),e2(),rho())-HS(underlying(),barrier(),2*mu())*M(f2(),-e4(),-rho())); 153 return result; 154 } else { 155 Real S1 = underlying()*std::exp((b-riskFreeRate())*residualTime()); 156 Real X1 = (strike()*std::exp(-riskFreeRate()*residualTime())); 157 Real HS1 = HS(underlying(),barrier(),2*(mu()+1)); 158 Real HS2 = HS(underlying(), barrier(), 2 * mu()); 159 result = S1; 160 result *= (M(-g1(),-e1(),rho())-HS1*M(-g3(),e3(),-rho())); 161 result -= X1*(M(-g2(), -e2(), rho()) - HS2*M(-g4(), e4(), -rho())); 162 result -= S1*(M(-d1(), -e1(), rho()) - HS1*M(-f1(), e3(), -rho())); 163 result += X1*(M(-d2(), -e2(), rho()) - HS2*M(-f2(), e4(), -rho())); 164 result += S1*(M(g1(), e1(), rho()) - HS1*M(g3(), -e3(), -rho())); 165 result -= X1*(M(g2(), e2(), rho()) - HS2*M(g4(), -e4(), -rho())); 166 return result; 167 } 168 } 169 170 // eta = -1: Up-and-In Call 171 // eta = 1: Down-and-In Call CIA(Integer eta) const172 Real AnalyticPartialTimeBarrierOptionEngine::CIA(Integer eta) const { 173 ext::shared_ptr<EuropeanExercise> exercise = 174 ext::dynamic_pointer_cast<EuropeanExercise>(arguments_.exercise); 175 176 ext::shared_ptr<PlainVanillaPayoff> payoff = 177 ext::dynamic_pointer_cast<PlainVanillaPayoff>(arguments_.payoff); 178 179 VanillaOption europeanOption(payoff, exercise); 180 181 europeanOption.setPricingEngine( 182 ext::make_shared<AnalyticEuropeanEngine>(process_)); 183 184 return europeanOption.NPV() - CA(eta); 185 } 186 CA(Integer eta) const187 Real AnalyticPartialTimeBarrierOptionEngine::CA(Integer eta) const { 188 //Partial-Time-Start- OUT Call Option calculation 189 Real b = riskFreeDiscount()-dividendYield(); 190 Real result; 191 result = underlying()*std::exp((b-riskFreeRate())*residualTime()); 192 result *= (M(d1(),eta*e1(),eta*rho())-HS(underlying(),barrier(),2*(mu()+1))*M(f1(),eta*e3(),eta*rho())); 193 result -= (strike()*std::exp(-riskFreeRate()*residualTime())*(M(d2(),eta*e2(),eta*rho())-HS(underlying(),barrier(),2*mu())*M(f2(),eta*e4(),eta*rho()))); 194 return result; 195 } 196 underlying() const197 Real AnalyticPartialTimeBarrierOptionEngine::underlying() const { 198 return process_->x0(); 199 } 200 strike() const201 Real AnalyticPartialTimeBarrierOptionEngine::strike() const { 202 ext::shared_ptr<PlainVanillaPayoff> payoff = 203 ext::dynamic_pointer_cast<PlainVanillaPayoff>(arguments_.payoff); 204 QL_REQUIRE(payoff, "non-plain payoff given"); 205 return payoff->strike(); 206 } 207 residualTime() const208 Time AnalyticPartialTimeBarrierOptionEngine::residualTime() const { 209 return process_->time(arguments_.exercise->lastDate()); 210 } 211 coverEventTime() const212 Time AnalyticPartialTimeBarrierOptionEngine::coverEventTime() const { 213 return process_->time(arguments_.coverEventDate); 214 } 215 volatility(Time t) const216 Volatility AnalyticPartialTimeBarrierOptionEngine::volatility(Time t) const { 217 return process_->blackVolatility()->blackVol(t, strike()); 218 } 219 stdDeviation() const220 Real AnalyticPartialTimeBarrierOptionEngine::stdDeviation() const { 221 Time T = residualTime(); 222 return volatility(T) * std::sqrt(T); 223 } 224 barrier() const225 Real AnalyticPartialTimeBarrierOptionEngine::barrier() const { 226 return arguments_.barrier; 227 } 228 rebate() const229 Real AnalyticPartialTimeBarrierOptionEngine::rebate() const { 230 return arguments_.rebate; 231 } 232 riskFreeRate() const233 Rate AnalyticPartialTimeBarrierOptionEngine::riskFreeRate() const { 234 return process_->riskFreeRate()->zeroRate(residualTime(), Continuous, 235 NoFrequency); 236 } 237 riskFreeDiscount() const238 DiscountFactor AnalyticPartialTimeBarrierOptionEngine::riskFreeDiscount() const { 239 return process_->riskFreeRate()->discount(residualTime()); 240 } 241 dividendYield() const242 Rate AnalyticPartialTimeBarrierOptionEngine::dividendYield() const { 243 return process_->dividendYield()->zeroRate(residualTime(), Continuous, 244 NoFrequency); 245 } 246 dividendDiscount() const247 DiscountFactor AnalyticPartialTimeBarrierOptionEngine::dividendDiscount() const { 248 return process_->dividendYield()->discount(residualTime()); 249 } 250 251 f1() const252 Real AnalyticPartialTimeBarrierOptionEngine::f1() const { 253 Real S = underlying(); 254 Real T = residualTime(); 255 Real sigma = volatility(T); 256 return (std::log(S / strike()) + 2 * std::log(barrier() / S) + ((riskFreeRate()-dividendYield()) + (std::pow(sigma, 2) / 2))*T) / (sigma*std::sqrt(T)); 257 } 258 f2() const259 Real AnalyticPartialTimeBarrierOptionEngine::f2() const { 260 Time T = residualTime(); 261 return f1() - volatility(T)*std::sqrt(T); 262 } 263 M(Real a,Real b,Real rho) const264 Real AnalyticPartialTimeBarrierOptionEngine::M(Real a,Real b,Real rho) const { 265 BivariateCumulativeNormalDistributionDr78 CmlNormDist(rho); 266 return CmlNormDist(a,b); 267 } 268 rho() const269 Real AnalyticPartialTimeBarrierOptionEngine::rho() const { 270 return std::sqrt(coverEventTime()/residualTime()); 271 } 272 mu() const273 Rate AnalyticPartialTimeBarrierOptionEngine::mu() const { 274 Volatility vol = volatility(coverEventTime()); 275 return ((riskFreeRate() - dividendYield()) - (vol * vol) / 2) / (vol * vol); 276 } 277 d1() const278 Real AnalyticPartialTimeBarrierOptionEngine::d1() const { 279 Real b = riskFreeRate()-dividendYield(); 280 Time T2 = residualTime(); 281 Volatility vol = volatility(T2); 282 return (std::log(underlying()/strike())+(b+vol*vol/2)*T2)/(std::sqrt(T2)*vol); 283 } 284 d2() const285 Real AnalyticPartialTimeBarrierOptionEngine::d2() const { 286 Time T2 = residualTime(); 287 Volatility vol = volatility(T2); 288 return d1() - vol*std::sqrt(T2); 289 } 290 e1() const291 Real AnalyticPartialTimeBarrierOptionEngine::e1() const { 292 Real b = riskFreeRate()-dividendYield(); 293 Time T1 = coverEventTime(); 294 Volatility vol = volatility(T1); 295 return (std::log(underlying()/barrier())+(b+vol*vol/2)*T1)/(std::sqrt(T1)*vol); 296 } 297 e2() const298 Real AnalyticPartialTimeBarrierOptionEngine::e2() const { 299 Time T1 = coverEventTime(); 300 Volatility vol = volatility(T1); 301 return e1() - vol*std::sqrt(T1); 302 } 303 e3() const304 Real AnalyticPartialTimeBarrierOptionEngine::e3() const { 305 Time T1 = coverEventTime(); 306 Real vol = volatility(T1); 307 return e1()+(2*std::log(barrier()/underlying()) /(vol*std::sqrt(T1))); 308 } 309 e4() const310 Real AnalyticPartialTimeBarrierOptionEngine::e4() const { 311 Time t = coverEventTime(); 312 return e3()-volatility(t)*std::sqrt(t); 313 } 314 g1() const315 Real AnalyticPartialTimeBarrierOptionEngine::g1() const { 316 Real b = riskFreeRate()-dividendYield(); 317 Time T2 = residualTime(); 318 Volatility vol = volatility(T2); 319 return (std::log(underlying()/barrier())+(b+vol*vol/2)*T2)/(std::sqrt(T2)*vol); 320 } 321 g2() const322 Real AnalyticPartialTimeBarrierOptionEngine::g2() const { 323 Time T2 = residualTime(); 324 Volatility vol = volatility(T2); 325 return g1() - vol*std::sqrt(T2); 326 } 327 g3() const328 Real AnalyticPartialTimeBarrierOptionEngine::g3() const { 329 Time T2 = residualTime(); 330 Real vol = volatility(T2); 331 return g1()+(2*std::log(barrier()/underlying()) /(vol*std::sqrt(T2))); 332 } 333 g4() const334 Real AnalyticPartialTimeBarrierOptionEngine::g4() const { 335 Time T2 = residualTime(); 336 Real vol = volatility(T2); 337 return g3()-vol*std::sqrt(T2); 338 } 339 HS(Real S,Real H,Real power) const340 Real AnalyticPartialTimeBarrierOptionEngine::HS(Real S, Real H, Real power) const { 341 return std::pow((H/S),power); 342 } 343 344 } 345 346