|
18 | 18 |
|
19 | 19 | /*! \defgroup findiff Finite-differences framework |
20 | 20 |
|
21 | | - <b>Warning: this section of the documentation is currently outdated.</b> |
22 | | - You will need to compare the information on this page with the present |
23 | | - code for working pricers, such as FdAmericanOption. |
24 | | - |
25 | 21 | This framework (corresponding to the ql/FiniteDifferences directory) |
26 | | - contains basic building blocks for the numerical solution of a generic |
27 | | - differential equation |
28 | | - \f[ |
29 | | - \frac{\partial f}{\partial t} = Lf |
30 | | - \f] |
31 | | - where \f$ L \f$ is a differential operator in ``space'', i.e., one which |
32 | | - does not contain partial derivatives in \f$ t \f$ but can otherwise |
33 | | - contain any derivative in any other variable of the problem. |
34 | | - |
35 | | - Writing the equation in the above form allows us to implement |
36 | | - separately the discretization of the differential operator \f$ L \f$ |
37 | | - and the time scheme used for the evolution of the solution. The |
38 | | - QuantLib::FiniteDifferenceModel class acts as a glue |
39 | | - for such two steps---which are outlined in the following sections---and |
40 | | - provides the interface of the resulting finite difference model for the |
41 | | - end user. Furthermore, it provides the possibility of checking and |
42 | | - operating on the solution array at each step---which is typically used |
43 | | - to apply an exercise condition for an option. This is also outlined in a |
44 | | - section below. |
45 | | - |
46 | | - |
47 | | - \section operators Differential operators |
48 | | - The discretization of the differential operator \f$ L \f$ depends on the |
49 | | - discretization chosen for the solution \f$ f \f$ of the given equation. |
50 | | - |
51 | | - Such choice is obvious in the 1-D case where the domain |
52 | | - \f$ [a,b] \f$ of the equation is discretized as a series of points |
53 | | - \f$ x_i, i=0 \dots N-1 \f$ (note that the index is zero based) where |
54 | | - \f$ x_i = a + hi \f$ and \f$ h = (b-a)/(N-1) \f$. In turn, the solution |
55 | | - \f$ f \f$ of the equation is discretized as an array |
56 | | - \f$ u_i, i=0 \dots N-1 \f$ whose elements are defined as |
57 | | - \f$ u_i = f(x_i) \f$. |
58 | | - The discretization of the differential operator follows by substituting |
59 | | - the derivatives with the corresponding incremental ratios defined in |
60 | | - terms of the \f$ f_i \f$. A number of basic operators are defined in the |
61 | | - framework which can be composed to form more complex operators, namely: |
62 | | - |
63 | | - the first derivative \f$ \partial/\partial x \f$ is discretized as the |
64 | | - operator \f$ D_+ \f$, defined as |
65 | | - \f[ D_{+} u_{i} = \frac{u_{i+1}-u_{i}}{h} \f] |
66 | | - and implemented in class QuantLib::DPlus; |
67 | | - the operator \f$ D_- \f$, defined as |
68 | | - \f[ D_{-} u_{i} = \frac{u_{i}-u_{i-1}}{h} \f] |
69 | | - and implemented in class QuantLib::DMinus; |
70 | | - and the operator \f$ D_0 \f$, defined as |
71 | | - \f[ D_{0} u_{i} = \frac{u_{i+1}-u_{i-1}}{2h} \f] |
72 | | - and implemented in class QuantLib::DZero. |
73 | | - The discretization error of the above operators is \f$ O(h) \f$ for |
74 | | - \f$ D_+ \f$ and \f$ D_- \f$ and \f$ O(h^2) \f$ for \f$ D_0 \f$; |
75 | | - |
76 | | - the second derivative \f$ \partial^2/\partial x^2 \f$ is |
77 | | - discretized as the operator \f$ D_+D_- \f$, defined as |
78 | | - \f[ D_{+}D_{-} u_{i} = \frac{u_{i+1}-2u_{i}+u_{i-1}}{h^2} \f] |
79 | | - and implemented in class QuantLib::DPlusDMinus. |
80 | | - Its discretization error is \f$ O(h^2) \f$. |
81 | | - |
82 | | - The boundary condition for the above operators is by default linear |
83 | | - extrapolation. Methods are currently provided for setting other kinds |
84 | | - of boundary conditions to a tridiagonal operator which these operators |
85 | | - inherit, namely, Dirichlet---i.e., constant value---and Neumann---i.e., |
86 | | - constant derivative---boundary conditions. |
87 | | - This might change in the future as boundary conditions could be |
88 | | - astracted and passed as an additional argument to the model. |
89 | | - |
90 | | - A programmer can also implement its own operator. |
91 | | - However, in order to fit into this framework it will have to implement |
92 | | - a required interface depending on the chosen evolver (see below). Also, |
93 | | - it is currently required to manage itself any boundary conditions. Again, |
94 | | - this could change in the future. |
95 | | - |
96 | | - On the other hand, there is no obvious choice in the 2-D case. |
97 | | - While it is immediate to discretize the domain into a series of points |
98 | | - \f$ (x_i,y_j) \f$ and the solution into a matrix |
99 | | - \f$ f_{ij} = f(x_i,y_j) \f$, there is a number of ways into which |
100 | | - the \f$ f_{ij} \f$ can be arranged into an array---each of them |
101 | | - determining a different discretization of the differential operators. |
102 | | - One of such ways was implemented in the LexicographicalView class, |
103 | | - while others will be implemented in the future. No 2-D operator is |
104 | | - currently implemented. |
105 | | - |
106 | | - |
107 | | - \section evolvers Time schemes |
108 | | - Once the differential operator \f$ L \f$ has been discretized, a number |
109 | | - of choices are available for discretizing the time derivative at the |
110 | | - left-hand side of the equation. |
111 | | - |
112 | | - In this framework, such choice is encapsulated in so-called evolvers |
113 | | - which, given \f$ L \f$ and the solution \f$ u^{(k)} \f$ at time |
114 | | - \f$ t_k \f$, yield the solution \f$ u^{(k-1)} \f$ at the previous |
115 | | - time step. |
116 | | - |
117 | | - A number of evolvers are currently provided in the library which |
118 | | - implement well-known schemes, namely, |
119 | | - |
120 | | - the forward Euler explicit scheme in which the equation is discretized as |
121 | | - \f[ \frac{u^{(k)}-u^{(k-1)}}{\Delta t} = Lu^{(k)} \f] |
122 | | - hence |
123 | | - \f[ u^{(k-1)} = \left( I - \Delta t L \right) u^{(k)} \f] |
124 | | - from which \f$ u^{(k-1)} \f$ can be obtained directly; |
125 | | - |
126 | | - the backward Euler implicit scheme in which the equation is discretized as |
127 | | - \f[ \frac{u^{(k)}-u^{(k-1)}}{\Delta t} = Lu^{(k-1)} \f] |
128 | | - hence |
129 | | - \f[ \left( I + \Delta t L \right) u^{(k-1)} = u^{(k)} \f] |
130 | | - from which \f$ u^{(k-1)} \f$ can be obtained by solving a linear system; |
131 | | - |
132 | | - the Crank-Nicolson scheme in which the equation is discretized as |
133 | | - \f[ \frac{u^{(k)}-u^{(k-1)}}{\Delta t} = |
134 | | - L \frac{u^{(k)}+u^{(k-1)}}{2} \f] |
135 | | - hence |
136 | | - \f[ \left( I + \frac{\Delta t}{2} L \right) u^{(k-1)} = |
137 | | - \left( I - \frac{\Delta t}{2} L \right) u^{(k)} \f] |
138 | | - from which \f$ u^{(k-1)} \f$ can be obtained by solving a linear system. |
139 | | - |
140 | | - Each of the above evolvers forces a set of interface requirements upon |
141 | | - the differential operator which are detailed in the documentation of the |
142 | | - corresponding class, namely, QuantLib::ExplicitEuler, |
143 | | - QuantLib::ImplicitEuler, and |
144 | | - QuantLib::CrankNicolson, respectively. |
145 | | - |
146 | | - A programmer could implement its own evolver, which does not need to |
147 | | - inherit from any base class. |
148 | | - |
149 | | - However, it must implement the following interface: |
150 | | - |
151 | | - \code |
152 | | - class Evolver { |
153 | | - public: |
154 | | - typedef ... arrayType; |
155 | | - typedef ... operatorType; |
156 | | - // constructors |
157 | | - Evolver(const operatorType& D); |
158 | | - // member functions |
159 | | - void step(arrayType& a, Time t) const; |
160 | | - void setStep(Time dt); |
161 | | - }; |
162 | | - \endcode |
163 | | - |
164 | | - Finally, we note again that the pricing of an option requires |
165 | | - the finite difference model to solve the corresponding equation |
166 | | - <em>backwards</em> in time. Therefore, given a discretization \f$ u \f$ |
167 | | - of the solution at a given time \f$ t \f$, the call |
168 | | - \code |
169 | | - evolver.step(u,t) |
170 | | - \endcode |
171 | | - must calculate the discrete solution at the <em>previous</em> time, |
172 | | - \f$ t-dt \f$. |
173 | | - |
174 | | - |
175 | | - \section conditions Step conditions |
176 | | - A finite difference model can be passed a step condition to be |
177 | | - applied at each step during the rollback of the solution (e.g. the early |
178 | | - exercise American condition). Such condition must be embodied in a class |
179 | | - derived from QuantLib::StepCondition and must implement |
180 | | - the interface of the latter, namely, |
181 | | - \code |
182 | | - class MyCondition : public StepCondition<arrayType> { |
183 | | - public: |
184 | | - void applyTo(arrayType& a, Time t) const; |
185 | | - }; |
186 | | - \endcode |
187 | | - |
188 | | - |
189 | | - \section fdexample An example of finite difference model |
190 | | - The Black-Scholes equation can be written in the above form as |
191 | | - \f[ \frac{\partial f}{\partial t} = |
192 | | - - \frac{\sigma^2}{2} \frac{\partial^2 f}{\partial x^2} |
193 | | - - \nu \frac{\partial f}{\partial x} |
194 | | - + r f. \f] |
195 | | - It can be seen that the operator \f$ L_{BS} \f$ is |
196 | | - \f[ L_{BS} = - \frac{\sigma^2}{2} \frac{\partial^2}{\partial x^2} |
197 | | - - \nu \frac{\partial}{\partial x} |
198 | | - + r I \f] |
199 | | - and can be built from the basic operators provided in the library as |
200 | | - \f[ L_{BS} = - \frac{\sigma^2}{2} D_{+}D_{-} |
201 | | - - \nu D_{0} + r I. \f] |
202 | | - |
203 | | - Its implementation closely reflects the above decomposition and |
204 | | - can be written as |
205 | | - \code |
206 | | - class BlackScholesOperator : public TridiagonalOperator { |
207 | | - public: |
208 | | - BlackScholesOperator( |
209 | | - double sigma, double nu, // parameters of the |
210 | | - Rate r, // Black-Scholes equation; |
211 | | - unsigned int points, // number of discretized points; |
212 | | - double h) // grid spacing. |
213 | | - : TridiagonalOperator( |
214 | | - // build the operator by adding basic ones |
215 | | - - (sigma*sigma/2.0) * DPlusDMinus(points,h) |
216 | | - - nu * DZero(points,h) |
217 | | - + r * TridiagonalOperator::identity(points) |
218 | | - ) {} |
219 | | - }; |
220 | | - \endcode |
221 | | - taking as inputs the relevant parameters of the equation |
222 | | - (\f$ \sigma \f$, \f$ \nu \f$ and \f$ r \f$) as well as model parameters |
223 | | - such as the number \f$ N \f$ of grid points and their spacing \f$ h \f$. |
224 | | - |
225 | | - As simple example cases, we will use the above operator to price both an |
226 | | - European and an American option. The parameters of the two options will |
227 | | - be the same, namely, they will be both call options with underlying price |
228 | | - \f$ u = 100 \f$, strike \f$ s = 95 \f$, residual time \f$ T = 1 \f$ year, |
229 | | - dividend yield \f$ q = 3\% \f$ and volatility \f$ \sigma = 10\% \f$. The |
230 | | - risk-free rate will be \f$ r = 5\% \f$. Such parameters are expressed |
231 | | - using %QuantLib types as |
232 | | - \code |
233 | | - Option::Type type = Option::Call; |
234 | | - double underlying = 100.0, strike = 95.0; |
235 | | - Time residualTime = 1.0; |
236 | | - Rate dividendYield = 0.03, riskFreeRate = 0.05; |
237 | | - double volatility = 0.10; |
238 | | - \endcode |
239 | | - |
240 | | - The grid upon which the model will act will be a logarithmic grid of |
241 | | - underlying prices, i.e., \f$ f \f$ will be defined in a range |
242 | | - \f$ [ \ln u_{min}, \ln u_{max}] \f$ discretized as an array |
243 | | - \f$ x_i, i = 0 \dots N-1 \f$ with \f$ x_i = \ln u_{min} + ih \f$ and |
244 | | - \f$ h = (\ln u_{max} - \ln u_{min})/(N-1) \f$. |
245 | | - Such a grid and the corresponding vector of actual prices can be built |
246 | | - as shown in the code below. The domain of the model will be defined |
247 | | - as \f$ [ \ln u - \Delta, \ln u + \Delta ] \f$ where |
248 | | - \f$ \Delta = 4 \sigma \sqrt{T} \f$. A number of grid points |
249 | | - \f$ N = 101 \f$ will be used. |
250 | | - \code |
251 | | - unsigned int gridPoints = 101; |
252 | | - Array grid(gridPoints), prices(gridPoints); |
253 | | - double x0 = QL_LOG(underlying); |
254 | | - double Delta = 4.0*volatility*QL_SQRT(residualTime); |
255 | | - double xMin = x0 - Delta, xMax = x0 + Delta; |
256 | | - double h = (xMax-xMin)/(gridPoints-1); |
257 | | - for (unsigned int i=0; i<gridPoints; i++) { |
258 | | - grid[i] = xMin + i*h; |
259 | | - prices[i] = QL_EXP(grid[i]); |
260 | | - } |
261 | | - \endcode |
262 | | - |
263 | | - The initial condition is determined by the values of the option at |
264 | | - maturity, i.e., either the difference between underlying price and |
265 | | - strike if such difference is positive, or 0 if that is not the case |
266 | | - (the above will have to be suitably modified for a put option or |
267 | | - a straddle.) |
268 | | - Such ``initial'' condition will be rolled back in time by our model. |
269 | | - \code |
270 | | - Array exercisingValue(gridPoints); |
271 | | - for (unsigned int i=0; i<gridPoints; i++) |
272 | | - exercisingValue[i] = QL_MAX(prices[i]-strike,0.0); |
273 | | - \endcode |
274 | | - |
275 | | - Now the differential operator can be initialized. Also, Neumann |
276 | | - initial conditions are set which correspond to the initial value |
277 | | - of the derivatives at the boundaries (see the BoundaryCondition |
278 | | - class documentation for details). |
279 | | - \code |
280 | | - double nu = riskFreeRate - dividendYield - volatility*volatility/2.0; |
281 | | - TridiagonalOperator L = BlackScholesOperator(volatility, nu, |
282 | | - riskFreeRate, gridPoints, h); |
283 | | - L.setLowerBC(BoundaryCondition(BoundaryCondition::Neumann, |
284 | | - exercisingValue[1]-exercisingValue[0])); |
285 | | - L.setUpperBC(BoundaryCondition(BoundaryCondition::Neumann, |
286 | | - exercisingValue[gridPoints_-1]-exercisingValue[gridPoints_-2])); |
287 | | - \endcode |
288 | | - |
289 | | - We are now already set for the pricing of the European option. |
290 | | - Also, the exercise condition is the only thing still to be defined |
291 | | - for the American option to be priced. Such condition is equivalent to |
292 | | - the statement that at each time step, the value of the option is the |
293 | | - maximum between the profit realized in exercising the option (which we |
294 | | - already calculated and stored in <tt>exercisingValue</tt>) and the value |
295 | | - of the option should we keep it (which corresponds to the solution rolled |
296 | | - back to the current time step). This logic can be implemented as: |
297 | | - \code |
298 | | - class ExerciseCondition : public StepCondition<Array> { |
299 | | - public: |
300 | | - ExerciseCondition(const Array& exercisingValue) |
301 | | - : exercisingValue_(exercisingValue) {} |
302 | | - void applyTo(Array& a, Time) const { |
303 | | - for (unsigned int i = 0; i < a.size(); i++) |
304 | | - a[i] = QL_MAX(a[i], exercisingValue_[i]); |
305 | | - } |
306 | | - private: |
307 | | - Array exercisingValue_; |
308 | | - }; |
309 | | - \endcode |
310 | | - |
311 | | - Everything is now ready. The model can be created gluing the piece |
312 | | - together by means of the QuantLib::FiniteDifferenceModel |
313 | | - class. The current value of the option is calculated by rolling back the |
314 | | - solution to the current time, i.e., \f$ t = 0 \f$, and by taking the value |
315 | | - corresponding at the current underlying price---which by construction |
316 | | - corresponds to the central value provided that the number of grid |
317 | | - points is odd. |
318 | | - \code |
319 | | - unsigned int timeSteps = 365; |
320 | | - |
321 | | - // build the model - Crank-Nicolson scheme chosen |
322 | | - FiniteDifferenceModel<CrankNicolson<TridiagonalOperator> > model(L); |
323 | | - |
324 | | - // European option |
325 | | - Array f = exercisingValue; // initial condition |
326 | | - model.rollback(f, residualTime, 0.0, timeSteps); |
327 | | - double europeanValue = valueAtCenter(f); |
| 22 | + contains basic building blocks for the numerical solution of partial |
| 23 | + differential equations by means of finite-difference methods. |
328 | 24 |
|
329 | | - // American option |
330 | | - f = exercisingValue; // reset |
331 | | - Handle<StepCondition<Array> > condition( |
332 | | - new ExerciseCondition(exercisingValue)); |
333 | | - model.rollback(f, residualTime, 0.0, timeSteps, condition); |
334 | | - double americanValue = valueAtCenter(f); |
335 | | - \endcode |
336 | 25 | */ |
0 commit comments