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/continuousarithmeticasianlevyengine.hpp>
21 #include <ql/pricingengines/blackcalculator.hpp>
22 #include <ql/math/distributions/normaldistribution.hpp>
23 #include <ql/exercise.hpp>
24 
25 using namespace std;
26 
27 namespace QuantLib {
28 
ContinuousArithmeticAsianLevyEngine(const ext::shared_ptr<GeneralizedBlackScholesProcess> & process,const Handle<Quote> & currentAverage,Date startDate)29     ContinuousArithmeticAsianLevyEngine::ContinuousArithmeticAsianLevyEngine(
30             const ext::shared_ptr<GeneralizedBlackScholesProcess>& process,
31             const Handle<Quote>& currentAverage,
32             Date startDate)
33     : process_(process), currentAverage_(currentAverage),
34       startDate_(startDate) {
35         registerWith(process_);
36         registerWith(currentAverage_);
37     }
38 
calculate() const39     void ContinuousArithmeticAsianLevyEngine::calculate() const {
40         QL_REQUIRE(arguments_.averageType == Average::Arithmetic,
41                    "not an Arithmetic average option");
42         QL_REQUIRE(arguments_.exercise->type() == Exercise::European,
43                    "not an European Option");
44         QL_REQUIRE(startDate_ <= process_->riskFreeRate()->referenceDate(),
45                    "startDate must be earlier than or equal to reference date");
46 
47         DayCounter rfdc  = process_->riskFreeRate()->dayCounter();
48         DayCounter divdc = process_->dividendYield()->dayCounter();
49         DayCounter voldc = process_->blackVolatility()->dayCounter();
50         Real spot = process_->stateVariable()->value();
51 
52         // payoff
53         ext::shared_ptr<StrikedTypePayoff> payoff =
54             ext::dynamic_pointer_cast<StrikedTypePayoff>(arguments_.payoff);
55         QL_REQUIRE(payoff, "non-plain payoff given");
56 
57         // original time to maturity
58         Date maturity = arguments_.exercise->lastDate();
59         Time T = rfdc.yearFraction(startDate_,
60                                    arguments_.exercise->lastDate());
61         // remaining time to maturity
62         Time T2 = rfdc.yearFraction(process_->riskFreeRate()->referenceDate(),
63                                     arguments_.exercise->lastDate());
64 
65         Real strike = payoff->strike();
66 
67         Volatility volatility =
68             process_->blackVolatility()->blackVol(maturity, strike);
69 
70         CumulativeNormalDistribution N;
71 
72         Rate riskFreeRate = process_->riskFreeRate()->
73             zeroRate(maturity, rfdc, Continuous, NoFrequency);
74         Rate dividendYield = process_->dividendYield()->
75             zeroRate(maturity, divdc, Continuous, NoFrequency);
76         Real b = riskFreeRate - dividendYield;
77 
78         Real Se = (std::fabs(b) > 1000*QL_EPSILON)
79             ? (spot/(T*b))*(exp((b-riskFreeRate)*T2)-exp(-riskFreeRate*T2))
80             : spot*T2/T * std::exp(-riskFreeRate*T2);
81 
82         Real X;
83         if (T2 < T) {
84             QL_REQUIRE(!currentAverage_.empty() && currentAverage_->isValid(),
85                        "current average required");
86             X = strike - ((T-T2)/T)*currentAverage_->value();
87         } else {
88             X = strike;
89         }
90 
91         Real m = (std::fabs(b) > 1000*QL_EPSILON) ? ((exp(b*T2)-1)/b) : T2;
92 
93         Real M = (2*spot*spot/(b+volatility*volatility)) *
94             (((exp((2*b+volatility*volatility)*T2)-1)
95               / (2*b+volatility*volatility))-m);
96 
97         Real D = M/(T*T);
98 
99         Real V = log(D)-2*(riskFreeRate*T2+log(Se));
100 
101         Real d1 = (1/sqrt(V))*((log(D)/2)-log(X));
102         Real d2 = d1-sqrt(V);
103 
104         if(payoff->optionType()==Option::Call)
105             results_.value = Se*N(d1) - X*exp(-riskFreeRate*T2)*N(d2);
106         else
107             results_.value = Se*N(d1) - X*exp(-riskFreeRate*T2)*N(d2)
108                              - Se + X*exp(-riskFreeRate*T2);
109     }
110 
111 }
112