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