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