1 /* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ 2 3 /* 4 Copyright (C) 2014 Bernd Lewerenz 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/continuousarithmeticasianvecerengine.hpp> 21 #include <ql/pricingengines/blackcalculator.hpp> 22 #include <ql/math/distributions/normaldistribution.hpp> 23 #include <ql/exercise.hpp> 24 #include <ql/math/rounding.hpp> 25 #include <ql/methods/finitedifferences/tridiagonaloperator.hpp> 26 #include <ql/methods/finitedifferences/dminus.hpp> 27 #include <ql/methods/finitedifferences/dplus.hpp> 28 #include <ql/methods/finitedifferences/dplusdminus.hpp> 29 #include <ql/instruments/vanillaoption.hpp> 30 #include <ql/pricingengines/vanilla/analyticeuropeanengine.hpp> 31 32 namespace QuantLib { 33 ContinuousArithmeticAsianVecerEngine(const ext::shared_ptr<GeneralizedBlackScholesProcess> & process,const Handle<Quote> & currentAverage,Date startDate,Size timeSteps,Size assetSteps,Real z_min,Real z_max)34 ContinuousArithmeticAsianVecerEngine::ContinuousArithmeticAsianVecerEngine( 35 const ext::shared_ptr<GeneralizedBlackScholesProcess>& process, 36 const Handle<Quote>& currentAverage, 37 Date startDate, 38 Size timeSteps, 39 Size assetSteps, 40 Real z_min, 41 Real z_max ) 42 : process_(process), currentAverage_(currentAverage), 43 startDate_(startDate),z_min_(z_min),z_max_(z_max), 44 timeSteps_(timeSteps),assetSteps_(assetSteps){ 45 registerWith(process_); 46 registerWith(currentAverage_); 47 } 48 calculate() const49 void ContinuousArithmeticAsianVecerEngine::calculate() const { 50 Real expectedAverage; 51 52 QL_REQUIRE(arguments_.averageType == Average::Arithmetic, 53 "not an Arithmetic average option"); 54 QL_REQUIRE(arguments_.exercise->type() == Exercise::European, 55 "not an European Option"); 56 57 DayCounter rfdc = process_->riskFreeRate()->dayCounter(); 58 DayCounter divdc = process_->dividendYield()->dayCounter(); 59 DayCounter voldc = process_->blackVolatility()->dayCounter(); 60 Real S_0 = process_->stateVariable()->value(); 61 62 // payoff 63 ext::shared_ptr<StrikedTypePayoff> payoff = 64 ext::dynamic_pointer_cast<StrikedTypePayoff>(arguments_.payoff); 65 QL_REQUIRE(payoff, "non-plain payoff given"); 66 67 // original time to maturity 68 Date maturity = arguments_.exercise->lastDate(); 69 70 Real X = payoff->strike(); 71 QL_REQUIRE(z_min_<=0 && z_max_>=0, 72 "strike (0 for vecer fixed strike asian) not on Grid"); 73 74 Volatility sigma = 75 process_->blackVolatility()->blackVol(maturity, X); 76 77 Rate r = process_->riskFreeRate()-> 78 zeroRate(maturity, rfdc, Continuous, NoFrequency); 79 Rate q = process_->dividendYield()-> 80 zeroRate(maturity, divdc, Continuous, NoFrequency); 81 82 Date today(Settings::instance().evaluationDate()); 83 84 QL_REQUIRE(startDate_>=today, 85 "Seasoned Asian not yet implemented"); 86 87 // Expiry in Years 88 Time T = rfdc.yearFraction(today, 89 arguments_.exercise->lastDate()); 90 Time T1 = rfdc.yearFraction(today, 91 startDate_ ); // Average Begin 92 Time T2 = T; // Average End (In this version only Maturity...) 93 94 if ((T2 - T1) < 0.001) { 95 // its a vanilla option. Use vanilla engine 96 VanillaOption europeanOption(payoff, arguments_.exercise); 97 europeanOption.setPricingEngine( 98 ext::make_shared<AnalyticEuropeanEngine>(process_)); 99 results_.value = europeanOption.NPV(); 100 101 } else { 102 Real Theta = 0.5; // Mixed Scheme: 0.5 = Crank Nicolson 103 Real Z_0 = cont_strategy(0,T1,T2,q,r) - std::exp(-r*T) * X /S_0; 104 105 QL_REQUIRE(Z_0>=z_min_ && Z_0<=z_max_, 106 "spot not on grid"); 107 108 Real h = (z_max_ - z_min_) / assetSteps_; // Space step size 109 Real k = T / timeSteps_; // Time Step size 110 111 Real sigma2 = sigma * sigma, vecerTerm; 112 113 Array SVec(assetSteps_+1),u_initial(assetSteps_+1), 114 u(assetSteps_+1),rhs(assetSteps_+1); 115 116 for (Natural i= 0; i<= SVec.size()-1;i++) { 117 SVec[i] = z_min_ + i * h; // Value of Underlying on the grid 118 } 119 120 // Begin gamma construction 121 TridiagonalOperator gammaOp = DPlusDMinus(assetSteps_+1,h); 122 123 Array upperD = gammaOp.upperDiagonal(); 124 Array lowerD = gammaOp.lowerDiagonal(); 125 Array Dia = gammaOp.diagonal(); 126 127 // Construct Vecer operator 128 TridiagonalOperator explicit_part(gammaOp.size()); 129 TridiagonalOperator implicit_part(gammaOp.size()); 130 131 for (Natural i= 0; i<= SVec.size()-1;i++) { 132 u_initial[i] = std::max<Real>(SVec[i] , 0.0); // Call Payoff 133 } 134 135 u = u_initial; 136 137 // Start Time Loop 138 139 for (Natural j = 1; j<=timeSteps_;j++) { 140 if (Theta != 1.0) { // Explicit Part 141 for (Natural i = 1; i<= SVec.size()-2;i++) { 142 vecerTerm = SVec[i] - std::exp(-q * (T-(j-1)*k)) 143 * cont_strategy(T-(j-1)*k,T1,T2,q,r); 144 gammaOp.setMidRow(i, 145 0.5 * sigma2 * vecerTerm * vecerTerm * lowerD[i-1], 146 0.5 * sigma2 * vecerTerm * vecerTerm * Dia[i], 147 0.5 * sigma2 * vecerTerm * vecerTerm * upperD[i]); 148 } 149 explicit_part = TridiagonalOperator::identity(gammaOp.size()) + 150 (1 - Theta) * k * gammaOp; 151 explicit_part.setFirstRow(1.0,0.0); // Apply before applying 152 explicit_part.setLastRow(-1.0,1.0); // Neumann BC 153 154 u = explicit_part.applyTo(u); 155 156 // Apply after applying (Neumann BC) 157 u[assetSteps_] = u[assetSteps_-1] + h; 158 u[0] = 0; 159 } // End Explicit Part 160 161 if (Theta != 0.0) { // Implicit Part 162 for (Natural i = 1; i<= SVec.size()-2;i++) { 163 vecerTerm = SVec[i] - std::exp(-q * (T-j*k)) * 164 cont_strategy(T-j*k,T1,T2,q,r); 165 gammaOp.setMidRow(i, 166 0.5 * sigma2 * vecerTerm * vecerTerm * lowerD[i-1], 167 0.5 * sigma2 * vecerTerm * vecerTerm * Dia[i], 168 0.5 * sigma2 * vecerTerm * vecerTerm * upperD[i]); 169 } 170 171 implicit_part = TridiagonalOperator::identity(gammaOp.size()) - 172 Theta * k * gammaOp; 173 174 // Apply before solving 175 implicit_part.setFirstRow(1.0,0.0); 176 implicit_part.setLastRow(-1.0,1.0); 177 rhs = u; 178 rhs[0] = 0; // Lower BC 179 rhs[assetSteps_] = h; // Upper BC (Neumann) Delta=1 180 u = implicit_part.solveFor(rhs); 181 } // End implicit Part 182 } // End Time Loop 183 184 DownRounding Rounding(0); 185 Integer lowerI = Integer(Rounding( (Z_0-z_min_)/h)); 186 // Interpolate solution 187 Real pv; 188 189 pv = u[lowerI] + (u[lowerI+1] - u[lowerI]) * (Z_0 - SVec[lowerI])/h; 190 results_.value = S_0 * pv; 191 192 if (payoff->optionType()==Option::Put) { 193 // Apply Call Put Parity for Asians 194 if (r == q) { 195 expectedAverage = S_0; 196 } else { 197 expectedAverage = 198 S_0 * (std::exp( (r-q) * T2) - 199 std::exp( (r-q) * T1)) / ((r-q) * (T2-T1)); 200 } 201 202 Real asianForward = std::exp(-r * T2) * (expectedAverage - X); 203 results_.value = results_.value - asianForward; 204 } 205 } 206 } 207 208 // Replication of Average by holding this amount in Assets cont_strategy(Time t,Time T1,Time T2,Real v,Real r) const209 Real ContinuousArithmeticAsianVecerEngine::cont_strategy(Time t, 210 Time T1, 211 Time T2, 212 Real v, 213 Real r) const { 214 Real const eps= 0.00001; 215 216 QL_REQUIRE(T1 <= T2, "Average Start must be before Average End"); 217 if (std::fabs(t-T2) < eps) { 218 return 0.0; 219 } else { 220 if (t<T1) { 221 if (std::fabs(r-v) >= eps) { 222 return (std::exp(v * (t-T2)) * 223 (1 - std::exp((v-r) * (T2-T1) )) / 224 (( r - v) * (T2 - T1) )); 225 } else { 226 return std::exp(v*(t-T2)); 227 } // end else v-r ==0 228 } else { // t<T1 229 if (std::fabs(r-v) >= eps) { 230 return std::exp(v * (t-T2)) * 231 (1 - std::exp( (v - r) * (T2-t) )) / 232 (( r - v) * (T2 - T1) ); 233 } else { 234 return std::exp(v * (t-T2)) * (T2 - t) / (T2 - T1); 235 } 236 } 237 } 238 } 239 240 } 241