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