Finite-differences framework
Detailed Description
Warning: this section of the documentation is currently outdated. You will need to compare the information on this page with the present code for working pricers, such as FdAmericanOption.This framework (corresponding to the ql/FiniteDifferences directory) contains basic building blocks for the numerical solution of a generic differential equation
where is a differential operator in ``space'', i.e., one which does not contain partial derivatives in
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 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

Such choice is obvious in the 1-D case where the domain of the equation is discretized as a series of points
(note that the index is zero based) where
and
. In turn, the solution
of the equation is discretized as an array
whose elements are defined as
. The discretization of the differential operator follows by substituting the derivatives with the corresponding incremental ratios defined in terms of the
. A number of basic operators are defined in the framework which can be composed to form more complex operators, namely:
the first derivative is discretized as the operator
, defined as
and implemented in class QuantLib::DPlus; the operator , defined as
and implemented in class QuantLib::DMinus; and the operator , defined as
and implemented in class QuantLib::DZero. The discretization error of the above operators is for
and
and
for
;
the second derivative is discretized as the operator
, defined as
and implemented in class QuantLib::DPlusDMinus. Its discretization error is .
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 and the solution into a matrix
, there is a number of ways into which the
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
In this framework, such choice is encapsulated in so-called evolvers which, given and the solution
at time
, yield the solution
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
hence
from which can be obtained directly;
the backward Euler implicit scheme in which the equation is discretized as
hence
from which can be obtained by solving a linear system;
the Crank-Nicolson scheme in which the equation is discretized as
hence
from which 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 of the solution at a given time
, the call
evolver.step(u,t)

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
It can be seen that the operator is
and can be built from the basic operators provided in the library as
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) ) {} };





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 , strike
, residual time
year, dividend yield
and volatility
. The risk-free rate will be
. 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., will be defined in a range
discretized as an array
with
and
. 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
where
. A number of grid points
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., , 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 |
![]() | |
class | DPlus |
![]() | |
class | DPlusDMinus |
![]() | |
class | DZero |
![]() | |
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 | OperatorFactory |
Black-Scholes-Merton differential operator. More... | |
class | StepConditionSet |
Parallel evolver for multiple arrays. More... | |
class | StepCondition |
condition to be applied at every time step More... | |
class | NullCondition |
null step condition More... | |
class | TridiagonalOperator |
Base implementation for tridiagonal operator. More... | |
Typedefs | |
typedef PdeOperator< PdeBSM > | QuantLib::BSMTermOperator |
Black-Scholes-Merton differential operator. | |
typedef PdeOperator< PdeShortRate > | QuantLib::OneFactorOperator |
Interest-rate single factor model differential operator. |
Typedef Documentation
|
Black-Scholes-Merton differential operator.
|