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