1 /* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ 2 3 /* 4 Copyright (C) 2011 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/analyticwriterextensibleoptionengine.hpp> 21 #include <ql/math/distributions/bivariatenormaldistribution.hpp> 22 #include <ql/pricingengines/blackformula.hpp> 23 24 using namespace std; 25 26 namespace QuantLib { 27 AnalyticWriterExtensibleOptionEngine(const ext::shared_ptr<GeneralizedBlackScholesProcess> & process)28 AnalyticWriterExtensibleOptionEngine::AnalyticWriterExtensibleOptionEngine( 29 const ext::shared_ptr<GeneralizedBlackScholesProcess>& process) 30 : process_(process) { 31 registerWith(process_); 32 } 33 calculate() const34 void AnalyticWriterExtensibleOptionEngine::calculate() const { 35 // We take all the arguments: 36 37 ext::shared_ptr<PlainVanillaPayoff> payoff1 = 38 ext::dynamic_pointer_cast<PlainVanillaPayoff>(arguments_.payoff); 39 QL_REQUIRE(payoff1, "not a plain vanilla payoff"); 40 41 ext::shared_ptr<PlainVanillaPayoff> payoff2 = 42 ext::dynamic_pointer_cast<PlainVanillaPayoff>(arguments_.payoff2); 43 QL_REQUIRE(payoff2, "not a plain vanilla payoff"); 44 45 ext::shared_ptr<Exercise> exercise1 = arguments_.exercise; 46 47 ext::shared_ptr<Exercise> exercise2 = arguments_.exercise2; 48 49 50 // We create and apply the calculate process: 51 52 Option::Type type = payoff1->optionType(); 53 54 // STEP 1: 55 56 // S = spot 57 Real spot = process_->stateVariable()->value(); 58 59 // For the B&S formulae: 60 DayCounter dividendDC = process_->dividendYield()->dayCounter(); 61 Rate dividend = process_->dividendYield()->zeroRate( 62 exercise1->lastDate(), dividendDC, Continuous, NoFrequency); 63 64 DayCounter riskFreeDC = process_->riskFreeRate()->dayCounter(); 65 Rate riskFree = process_->riskFreeRate()->zeroRate( 66 exercise1->lastDate(), riskFreeDC, Continuous, NoFrequency); 67 68 // The time to maturity: 69 Time t1 = riskFreeDC.yearFraction( 70 process_->riskFreeRate()->referenceDate(), 71 arguments_.exercise->lastDate()); 72 Time t2 = riskFreeDC.yearFraction( 73 process_->riskFreeRate()->referenceDate(), 74 arguments_.exercise2->lastDate()); 75 76 // b = r-q: 77 Real b = riskFree - dividend; 78 79 Real forwardPrice = spot * std::exp(b*t1); 80 81 Volatility volatility = process_->blackVolatility()->blackVol( 82 exercise1->lastDate(), payoff1->strike()); 83 84 Real stdDev = volatility*std::sqrt(t1); 85 86 Real discount = std::exp(-riskFree*t1); 87 88 // Call the B&S method: 89 Real black = blackFormula(type, payoff1->strike(), 90 forwardPrice, stdDev, discount); 91 92 // STEP 2: 93 94 // Standard bivariate normal distribution: 95 Real ro = std::sqrt(t1/t2); 96 Real z1 = (std::log(spot/payoff2->strike()) + 97 (b+std::pow(volatility, 2)/2)*t2)/(volatility*std::sqrt(t2)); 98 Real z2 = (std::log(spot/payoff1->strike()) + 99 (b+std::pow(volatility, 2)/2)*t1)/(volatility*std::sqrt(t1)); 100 101 // Call the bivariate method: 102 BivariateCumulativeNormalDistributionWe04DP biv(-ro); 103 104 105 // STEP 3: 106 107 Real bivariate1, bivariate2, result; 108 109 // Final computing: 110 if (type == Option::Call) { 111 // Call case: 112 bivariate1 = biv(z1, -z2); 113 bivariate2 = biv(z1-volatility*std::sqrt(t2), 114 -z2+volatility*std::sqrt(t1)); 115 result = black + spot*std::exp((b-riskFree)*t2)*bivariate1 116 - payoff2->strike()*std::exp((-riskFree)*t2)*bivariate2; 117 } else { 118 // Put case: 119 bivariate1 = biv(-z1, z2); 120 bivariate2 = biv(-z1+volatility*std::sqrt(t2), 121 z2-volatility*std::sqrt(t1)); 122 result = black - spot*std::exp((b-riskFree)*t2)*bivariate1 123 + payoff2->strike()*std::exp((-riskFree)*t2)*bivariate2; 124 } 125 126 // Save the result: 127 results_.value = result; 128 } 129 130 } 131