Finite-differences framework


Detailed Description

This framework (corresponding to the ql/FiniteDifferences directory) contains basic building blocks for the numerical solution of a generic differential equation

\[ \frac{\partial f}{\partial t} = Lf \]

where $ L $ is a differential operator in ``space'', i.e., one which does not contain partial derivatives in $ t $ but can otherwise contain any derivative in any other variable of the problem.

Writing the equation in the above form allows us to implement separately the discretization of the differential operator $ L $ and the time scheme used for the evolution of the solution. The QuantLib::FiniteDifferenceModel class acts as a glue for such two steps---which are outlined in the following sections---and provides the interface of the resulting finite difference model for the end user. Furthermore, it provides the possibility of checking and operating on the solution array at each step---which is typically used to apply an exercise condition for an option. This is also outlined in a section below.

Differential operators

The discretization of the differential operator $ L $ depends on the discretization chosen for the solution $ f $ of the given equation.

Such choice is obvious in the 1-D case where the domain $ [a,b] $ of the equation is discretized as a series of points $ x_i, i=0 \dots N-1 $ (note that the index is zero based) where $ x_i = a + hi $ and $ h = (b-a)/(N-1) $. In turn, the solution $ f $ of the equation is discretized as an array $ u_i, i=0 \dots N-1 $ whose elements are defined as $ u_i = f(x_i) $. The discretization of the differential operator follows by substituting the derivatives with the corresponding incremental ratios defined in terms of the $ f_i $. A number of basic operators are defined in the framework which can be composed to form more complex operators, namely:

the first derivative $ \partial/\partial x $ is discretized as the operator $ D_+ $, defined as

\[ D_{+} u_{i} = \frac{u_{i+1}-u_{i}}{h} \]

and implemented in class QuantLib::DPlus; the operator $ D_- $, defined as

\[ D_{-} u_{i} = \frac{u_{i}-u_{i-1}}{h} \]

and implemented in class QuantLib::DMinus; and the operator $ D_0 $, defined as

\[ D_{0} u_{i} = \frac{u_{i+1}-u_{i-1}}{2h} \]

and implemented in class QuantLib::DZero. The discretization error of the above operators is $ O(h) $ for $ D_+ $ and $ D_- $ and $ O(h^2) $ for $ D_0 $;

the second derivative $ \partial^2/\partial x^2 $ is discretized as the operator $ D_+D_- $, defined as

\[ D_{+}D_{-} u_{i} = \frac{u_{i+1}-2u_{i}+u_{i-1}}{h^2} \]

and implemented in class QuantLib::DPlusDMinus. Its discretization error is $ O(h^2) $.

The boundary condition for the above operators is by default linear extrapolation. Methods are currently provided for setting other kinds of boundary conditions to a tridiagonal operator which these operators inherit, namely, Dirichlet---i.e., constant value---and Neumann---i.e., constant derivative---boundary conditions. This might change in the future as boundary conditions could be astracted and passed as an additional argument to the model.

A programmer can also implement its own operator. However, in order to fit into this framework it will have to implement a required interface depending on the chosen evolver (see below). Also, it is currently required to manage itself any boundary conditions. Again, this could change in the future.

On the other hand, there is no obvious choice in the 2-D case. While it is immediate to discretize the domain into a series of points $ (x_i,y_j) $ and the solution into a matrix $ f_{ij} = f(x_i,y_j) $, there is a number of ways into which the $ f_{ij} $ can be arranged into an array---each of them determining a different discretization of the differential operators. One of such ways was implemented in the LexicographicalView class, while others will be implemented in the future. No 2-D operator is currently implemented.

Time schemes

Once the differential operator $ L $ has been discretized, a number of choices are available for discretizing the time derivative at the left-hand side of the equation.

In this framework, such choice is encapsulated in so-called evolvers which, given $ L $ and the solution $ u^{(k)} $ at time $ t_k $, yield the solution $ u^{(k-1)} $ at the previous time step.

A number of evolvers are currently provided in the library which implement well-known schemes, namely,

the forward Euler explicit scheme in which the equation is discretized as

\[ \frac{u^{(k)}-u^{(k-1)}}{\Delta t} = Lu^{(k)} \]

hence

\[ u^{(k-1)} = \left( I - \Delta t L \right) u^{(k)} \]

from which $ u^{(k-1)} $ can be obtained directly;

the backward Euler implicit scheme in which the equation is discretized as

\[ \frac{u^{(k)}-u^{(k-1)}}{\Delta t} = Lu^{(k-1)} \]

hence

\[ \left( I + \Delta t L \right) u^{(k-1)} = u^{(k)} \]

from which $ u^{(k-1)} $ can be obtained by solving a linear system;

the Crank-Nicolson scheme in which the equation is discretized as

\[ \frac{u^{(k)}-u^{(k-1)}}{\Delta t} = L \frac{u^{(k)}+u^{(k-1)}}{2} \]

hence

\[ \left( I + \frac{\Delta t}{2} L \right) u^{(k-1)} = \left( I - \frac{\Delta t}{2} L \right) u^{(k)} \]

from which $ u^{(k-1)} $ can be obtained by solving a linear system.

Each of the above evolvers forces a set of interface requirements upon the differential operator which are detailed in the documentation of the corresponding class, namely, QuantLib::ExplicitEuler, QuantLib::ImplicitEuler, and QuantLib::CrankNicolson, respectively.

A programmer could implement its own evolver, which does not need to inherit from any base class.

However, it must implement the following interface:

class Evolver { public: typedef ... arrayType; typedef ... operatorType; // constructors Evolver(const operatorType& D); // member functions void step(arrayType& a, Time t) const; void setStep(Time dt); };

Finally, we note again that the pricing of an option requires the finite difference model to solve the corresponding equation backwards in time. Therefore, given a discretization $ u $ of the solution at a given time $ t $, the call

evolver.step(u,t)
must calculate the discrete solution at the previous time, $ t-dt $.

Step conditions

A finite difference model can be passed a step condition to be applied at each step during the rollback of the solution (e.g. the early exercise American condition). Such condition must be embodied in a class derived from QuantLib::StepCondition and must implement the interface of the latter, namely,
class MyCondition : public StepCondition<arrayType> { public: void applyTo(arrayType& a, Time t) const; };

An example of finite difference model

The Black-Scholes equation can be written in the above form as

\[ \frac{\partial f}{\partial t} = - \frac{\sigma^2}{2} \frac{\partial^2 f}{\partial x^2} - \nu \frac{\partial f}{\partial x} + r f. \]

It can be seen that the operator $ L_{BS} $ is

\[ L_{BS} = - \frac{\sigma^2}{2} \frac{\partial^2}{\partial x^2} - \nu \frac{\partial}{\partial x} + r I \]

and can be built from the basic operators provided in the library as

\[ L_{BS} = - \frac{\sigma^2}{2} D_{+}D_{-} - \nu D_{0} + r I. \]

Its implementation closely reflects the above decomposition and can be written as

class BlackScholesOperator : public TridiagonalOperator { public: BlackScholesOperator( double sigma, double nu, // parameters of the Rate r, // Black-Scholes equation; unsigned int points, // number of discretized points; double h) // grid spacing. : TridiagonalOperator( // build the operator by adding basic ones - (sigma*sigma/2.0) * DPlusDMinus(points,h) - nu * DZero(points,h) + r * TridiagonalOperator::identity(points) ) {} };
taking as inputs the relevant parameters of the equation ($ \sigma $, $ \nu $ and $ r $) as well as model parameters such as the number $ N $ of grid points and their spacing $ h $.

As simple example cases, we will use the above operator to price both an European and an American option. The parameters of the two options will be the same, namely, they will be both call options with underlying price $ u = 100 $, strike $ s = 95 $, residual time $ T = 1 $ year, dividend yield $ q = 3\% $ and volatility $ \sigma = 10\% $. The risk-free rate will be $ r = 5\% $. Such parameters are expressed using QuantLib types as

Option::Type type = Option::Call; double underlying = 100.0, strike = 95.0; Time residualTime = 1.0; Rate dividendYield = 0.03, riskFreeRate = 0.05; double volatility = 0.10;

The grid upon which the model will act will be a logarithmic grid of underlying prices, i.e., $ f $ will be defined in a range $ [ \ln u_{min}, \ln u_{max}] $ discretized as an array $ x_i, i = 0 \dots N-1 $ with $ x_i = \ln u_{min} + ih $ and $ h = (\ln u_{max} - \ln u_{min})/(N-1) $. Such a grid and the corresponding vector of actual prices can be built as shown in the code below. The domain of the model will be defined as $ [ \ln u - \Delta, \ln u + \Delta ] $ where $ \Delta = 4 \sigma \sqrt{T} $. A number of grid points $ N = 101 $ will be used.

unsigned int gridPoints = 101; Array grid(gridPoints), prices(gridPoints); double x0 = QL_LOG(underlying); double Delta = 4.0*volatility*QL_SQRT(residualTime); double xMin = x0 - Delta, xMax = x0 + Delta; double h = (xMax-xMin)/(gridPoints-1); for (unsigned int i=0; i<gridPoints; i++) { grid[i] = xMin + i*h; prices[i] = QL_EXP(grid[i]); }

The initial condition is determined by the values of the option at maturity, i.e., either the difference between underlying price and strike if such difference is positive, or 0 if that is not the case (the above will have to be suitably modified for a put option or a straddle.) Such ``initial'' condition will be rolled back in time by our model.

Array exercisingValue(gridPoints); for (unsigned int i=0; i<gridPoints; i++) exercisingValue[i] = QL_MAX(prices[i]-strike,0.0);

Now the differential operator can be initialized. Also, Neumann initial conditions are set which correspond to the initial value of the derivatives at the boundaries (see the BoundaryCondition class documentation for details).

double nu = riskFreeRate - dividendYield - volatility*volatility/2.0; TridiagonalOperator L = BlackScholesOperator(volatility, nu, riskFreeRate, gridPoints, h); L.setLowerBC(BoundaryCondition(BoundaryCondition::Neumann, exercisingValue[1]-exercisingValue[0])); L.setUpperBC(BoundaryCondition(BoundaryCondition::Neumann, exercisingValue[gridPoints_-1]-exercisingValue[gridPoints_-2]));

We are now already set for the pricing of the European option. Also, the exercise condition is the only thing still to be defined for the American option to be priced. Such condition is equivalent to the statement that at each time step, the value of the option is the maximum between the profit realized in exercising the option (which we already calculated and stored in exercisingValue) and the value of the option should we keep it (which corresponds to the solution rolled back to the current time step). This logic can be implemented as:

class ExerciseCondition : public StepCondition<Array> { public: ExerciseCondition(const Array& exercisingValue) : exercisingValue_(exercisingValue) {} void applyTo(Array& a, Time) const { for (unsigned int i = 0; i < a.size(); i++) a[i] = QL_MAX(a[i], exercisingValue_[i]); } private: Array exercisingValue_; };

Everything is now ready. The model can be created gluing the piece together by means of the QuantLib::FiniteDifferenceModel class. The current value of the option is calculated by rolling back the solution to the current time, i.e., $ t = 0 $, and by taking the value corresponding at the current underlying price---which by construction corresponds to the central value provided that the number of grid points is odd.

unsigned int timeSteps = 365; // build the model - Crank-Nicolson scheme chosen FiniteDifferenceModel<CrankNicolson<TridiagonalOperator> > model(L); // European option Array f = exercisingValue; // initial condition model.rollback(f, residualTime, 0.0, timeSteps); double europeanValue = valueAtCenter(f); // American option f = exercisingValue; // reset Handle<StepCondition<Array> > condition( new ExerciseCondition(exercisingValue)); model.rollback(f, residualTime, 0.0, timeSteps, condition); double americanValue = valueAtCenter(f);


Classes

class  BoundaryCondition
 Abstract boundary condition class for finite difference problems. More...

class  NeumannBC
 Neumann boundary condition (i.e., constant derivative). More...

class  DirichletBC
 Neumann boundary condition (i.e., constant value). More...

class  BSMOperator
 Black-Scholes-Merton differential operator. More...

class  CrankNicolson
 Crank-Nicolson scheme for finite difference methods. More...

class  DMinus
 $ D_{-} $ matricial representation More...

class  DPlus
 $ D_{+} $ matricial representation More...

class  DPlusDMinus
 $ D_{+}D_{-} $ matricial representation More...

class  DZero
 $ D_{0} $ matricial representation More...

class  ExplicitEuler
 Forward Euler scheme for finite difference methods. More...

class  FiniteDifferenceModel
 Generic finite difference model. More...

class  ImplicitEuler
 Backward Euler scheme for finite difference methods. More...

class  MixedScheme
 Mixed (explicit/implicit) scheme for finite difference methods. More...

class  OneFactorOperator
 Interest-rate single factor model differential operator. More...

class  StepCondition
 condition to be applied at every time step More...

class  TridiagonalOperator
 Base implementation for tridiagonal operator. More...


QuantLib.org
QuantLib
Hosted by
SourceForge.net Logo
Documentation generated by
doxygen