QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
continuousarithmeticasianvecerengine.cpp
Go to the documentation of this file.
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/exercise.hpp>
24#include <ql/math/rounding.hpp>
31#include <utility>
32
33namespace QuantLib {
34
36 ext::shared_ptr<GeneralizedBlackScholesProcess> process,
37 Handle<Quote> currentAverage,
38 Date startDate,
39 Size timeSteps,
40 Size assetSteps,
41 Real z_min,
42 Real z_max)
43 : process_(std::move(process)), currentAverage_(std::move(currentAverage)),
44 startDate_(startDate), z_min_(z_min), z_max_(z_max), timeSteps_(timeSteps),
45 assetSteps_(assetSteps) {
48 }
49
51 Real expectedAverage;
52
54 "not an Arithmetic average option");
55 QL_REQUIRE(arguments_.exercise->type() == Exercise::European,
56 "not an European Option");
57
58 DayCounter rfdc = process_->riskFreeRate()->dayCounter();
59 DayCounter divdc = process_->dividendYield()->dayCounter();
60 DayCounter voldc = process_->blackVolatility()->dayCounter();
61 Real S_0 = process_->stateVariable()->value();
62
63 // payoff
64 ext::shared_ptr<StrikedTypePayoff> payoff =
65 ext::dynamic_pointer_cast<StrikedTypePayoff>(arguments_.payoff);
66 QL_REQUIRE(payoff, "non-plain payoff given");
67
68 // original time to maturity
69 Date maturity = arguments_.exercise->lastDate();
70
71 Real X = payoff->strike();
72 QL_REQUIRE(z_min_<=0 && z_max_>=0,
73 "strike (0 for vecer fixed strike asian) not on Grid");
74
76 process_->blackVolatility()->blackVol(maturity, X);
77
78 Rate r = process_->riskFreeRate()->
79 zeroRate(maturity, rfdc, Continuous, NoFrequency);
80 Rate q = process_->dividendYield()->
81 zeroRate(maturity, divdc, Continuous, NoFrequency);
82
83 Date today(Settings::instance().evaluationDate());
84
86 "Seasoned Asian not yet implemented");
87
88 // Expiry in Years
89 Time T = rfdc.yearFraction(today,
90 arguments_.exercise->lastDate());
91 Time T1 = rfdc.yearFraction(today,
92 startDate_ ); // Average Begin
93 Time T2 = T; // Average End (In this version only Maturity...)
94
95 if ((T2 - T1) < 0.001) {
96 // its a vanilla option. Use vanilla engine
97 VanillaOption europeanOption(payoff, arguments_.exercise);
98 europeanOption.setPricingEngine(
99 ext::make_shared<AnalyticEuropeanEngine>(process_));
100 results_.value = europeanOption.NPV();
101
102 } else {
103 Real Theta = 0.5; // Mixed Scheme: 0.5 = Crank Nicolson
104 Real Z_0 = cont_strategy(0,T1,T2,q,r) - std::exp(-r*T) * X /S_0;
105
106 QL_REQUIRE(Z_0>=z_min_ && Z_0<=z_max_,
107 "spot not on grid");
108
109 Real h = (z_max_ - z_min_) / assetSteps_; // Space step size
110 Real k = T / timeSteps_; // Time Step size
111
112 Real sigma2 = sigma * sigma, vecerTerm;
113
114 Array SVec(assetSteps_+1),u_initial(assetSteps_+1),
115 u(assetSteps_+1),rhs(assetSteps_+1);
116
117 for (Natural i= 0; i<= SVec.size()-1;i++) {
118 SVec[i] = z_min_ + i * h; // Value of Underlying on the grid
119 }
120
121 // Begin gamma construction
123
124 Array upperD = gammaOp.upperDiagonal();
125 Array lowerD = gammaOp.lowerDiagonal();
126 Array Dia = gammaOp.diagonal();
127
128 // Construct Vecer operator
129 TridiagonalOperator explicit_part(gammaOp.size());
130 TridiagonalOperator implicit_part(gammaOp.size());
131
132 for (Natural i= 0; i<= SVec.size()-1;i++) {
133 u_initial[i] = std::max<Real>(SVec[i] , 0.0); // Call Payoff
134 }
135
136 u = u_initial;
137
138 // Start Time Loop
139
140 for (Natural j = 1; j<=timeSteps_;j++) {
141 if (Theta != 1.0) { // Explicit Part
142 for (Natural i = 1; i<= SVec.size()-2;i++) {
143 vecerTerm = SVec[i] - std::exp(-q * (T-(j-1)*k))
144 * cont_strategy(T-(j-1)*k,T1,T2,q,r);
145 gammaOp.setMidRow(i,
146 0.5 * sigma2 * vecerTerm * vecerTerm * lowerD[i-1],
147 0.5 * sigma2 * vecerTerm * vecerTerm * Dia[i],
148 0.5 * sigma2 * vecerTerm * vecerTerm * upperD[i]);
149 }
150 explicit_part = TridiagonalOperator::identity(gammaOp.size()) +
151 (1 - Theta) * k * gammaOp;
152 explicit_part.setFirstRow(1.0,0.0); // Apply before applying
153 explicit_part.setLastRow(-1.0,1.0); // Neumann BC
154
155 u = explicit_part.applyTo(u);
156
157 // Apply after applying (Neumann BC)
158 u[assetSteps_] = u[assetSteps_-1] + h;
159 u[0] = 0;
160 } // End Explicit Part
161
162 if (Theta != 0.0) { // Implicit Part
163 for (Natural i = 1; i<= SVec.size()-2;i++) {
164 vecerTerm = SVec[i] - std::exp(-q * (T-j*k)) *
165 cont_strategy(T-j*k,T1,T2,q,r);
166 gammaOp.setMidRow(i,
167 0.5 * sigma2 * vecerTerm * vecerTerm * lowerD[i-1],
168 0.5 * sigma2 * vecerTerm * vecerTerm * Dia[i],
169 0.5 * sigma2 * vecerTerm * vecerTerm * upperD[i]);
170 }
171
172 implicit_part = TridiagonalOperator::identity(gammaOp.size()) -
173 Theta * k * gammaOp;
174
175 // Apply before solving
176 implicit_part.setFirstRow(1.0,0.0);
177 implicit_part.setLastRow(-1.0,1.0);
178 rhs = u;
179 rhs[0] = 0; // Lower BC
180 rhs[assetSteps_] = h; // Upper BC (Neumann) Delta=1
181 u = implicit_part.solveFor(rhs);
182 } // End implicit Part
183 } // End Time Loop
184
186 auto lowerI = Integer(Rounding((Z_0 - z_min_) / h));
187 // Interpolate solution
188 Real pv;
189
190 pv = u[lowerI] + (u[lowerI+1] - u[lowerI]) * (Z_0 - SVec[lowerI])/h;
191 results_.value = S_0 * pv;
192
193 if (payoff->optionType()==Option::Put) {
194 // Apply Call Put Parity for Asians
195 if (r == q) {
196 expectedAverage = S_0;
197 } else {
198 expectedAverage =
199 S_0 * (std::exp( (r-q) * T2) -
200 std::exp( (r-q) * T1)) / ((r-q) * (T2-T1));
201 }
202
203 Real asianForward = std::exp(-r * T2) * (expectedAverage - X);
204 results_.value = results_.value - asianForward;
205 }
206 }
207 }
208
209 // Replication of Average by holding this amount in Assets
211 Time T1,
212 Time T2,
213 Real v,
214 Real r) const {
215 Real const eps= 0.00001;
216
217 QL_REQUIRE(T1 <= T2, "Average Start must be before Average End");
218 if (std::fabs(t-T2) < eps) {
219 return 0.0;
220 } else {
221 if (t<T1) {
222 if (std::fabs(r-v) >= eps) {
223 return (std::exp(v * (t-T2)) *
224 (1 - std::exp((v-r) * (T2-T1) )) /
225 (( r - v) * (T2 - T1) ));
226 } else {
227 return std::exp(v*(t-T2));
228 } // end else v-r ==0
229 } else { // t<T1
230 if (std::fabs(r-v) >= eps) {
231 return std::exp(v * (t-T2)) *
232 (1 - std::exp( (v - r) * (T2-t) )) /
233 (( r - v) * (T2 - T1) );
234 } else {
235 return std::exp(v * (t-T2)) * (T2 - t) / (T2 - T1);
236 }
237 }
238 }
239 }
240
241}
Analytic European engine.
Black-formula calculator class.
1-D array used in linear algebra.
Definition: array.hpp:52
Real cont_strategy(Time t, Time T1, Time T2, Real v, Real r) const
ContinuousArithmeticAsianVecerEngine(ext::shared_ptr< GeneralizedBlackScholesProcess > process, Handle< Quote > currentAverage, Date startDate, Size timeSteps=100, Size assetSteps=100, Real z_min=-1.0, Real z_max=1.0)
ext::shared_ptr< GeneralizedBlackScholesProcess > process_
matricial representation
Definition: dplusdminus.hpp:43
Concrete date class.
Definition: date.hpp:125
day counter class
Definition: daycounter.hpp:44
Time yearFraction(const Date &, const Date &, const Date &refPeriodStart=Date(), const Date &refPeriodEnd=Date()) const
Returns the period between two dates as a fraction of year.
Definition: daycounter.hpp:128
Down-rounding.
Definition: rounding.hpp:98
Shared handle to an observable.
Definition: handle.hpp:41
Real NPV() const
returns the net present value of the instrument.
Definition: instrument.hpp:167
void setPricingEngine(const ext::shared_ptr< PricingEngine > &)
set the pricing engine to be used.
Definition: instrument.cpp:35
std::pair< iterator, bool > registerWith(const ext::shared_ptr< Observable > &)
Definition: observable.hpp:228
basic rounding class
Definition: rounding.hpp:35
static Settings & instance()
access to the unique instance
Definition: singleton.hpp:104
Base implementation for tridiagonal operator.
Array solveFor(const Array &rhs) const
solve linear system for a given right-hand side
Array applyTo(const Array &v) const
apply operator to a given array
const Array & lowerDiagonal() const
const Array & upperDiagonal() const
static TridiagonalOperator identity(Size size)
identity instance
void setMidRow(Size, Real, Real, Real)
Vanilla option (no discrete dividends, no barriers) on a single asset.
Vecer engine for continuous arithmetic Asian options.
const DefaultType & t
matricial representation
matricial representation
matricial representation
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
Definition: errors.hpp:117
Option exercise classes and payoff function.
@ NoFrequency
null frequency
Definition: frequency.hpp:37
Real Time
continuous quantity with 1-year units
Definition: types.hpp:62
QL_REAL Real
real number
Definition: types.hpp:50
unsigned QL_INTEGER Natural
positive integer
Definition: types.hpp:43
Real Volatility
volatility
Definition: types.hpp:78
QL_INTEGER Integer
integer number
Definition: types.hpp:35
Real Rate
interest rates
Definition: types.hpp:70
std::size_t Size
size of a container
Definition: types.hpp:58
Real sigma
ext::shared_ptr< QuantLib::Payoff > payoff
Definition: any.hpp:35
STL namespace.
normal, cumulative and inverse cumulative distributions
ext::shared_ptr< YieldTermStructure > q
ext::shared_ptr< YieldTermStructure > r
ext::shared_ptr< BlackVolTermStructure > v
Rounding implementation.
tridiagonal operator
Vanilla option on a single asset.