NAG Library Routine Document
D03NCF
1 Purpose
D03NCF solves the Black–Scholes equation for financial option pricing using a finite difference scheme.
2 Specification
SUBROUTINE D03NCF ( |
KOPT, X, MESH, NS, S, NT, T, TDPAR, R, Q, SIGMA, ALPHA, NTKEEP, F, THETA, DELTA, GAMMA, LAMBDA, RHO, LDF, WORK, IWORK, IFAIL) |
INTEGER |
KOPT, NS, NT, NTKEEP, LDF, IWORK(NS), IFAIL |
REAL (KIND=nag_wp) |
X, S(NS), T(NT), R(*), Q(*), SIGMA(*), ALPHA, F(LDF,NTKEEP), THETA(LDF,NTKEEP), DELTA(LDF,NTKEEP), GAMMA(LDF,NTKEEP), LAMBDA(LDF,NTKEEP), RHO(LDF,NTKEEP), WORK(4*NS) |
LOGICAL |
TDPAR(3) |
CHARACTER(1) |
MESH |
|
3 Description
D03NCF solves the Black–Scholes equation (see
Hull (1989) and
Wilmott et al. (1995))
for the value
of a European or American, put or call stock option, with exercise price
. In equation
(1) is time,
is the stock price,
is the risk free interest rate,
is the continuous dividend, and
is the stock volatility. According to the values in the array
TDPAR, the parameters
,
and
may each be either constant or functions of time. The routine also returns values of various Greeks.
D03NCF uses a finite difference method with a choice of time-stepping schemes. The method is explicit for
and implicit for nonzero values of
ALPHA. Second order time accuracy can be obtained by setting
. According to the value of the parameter
MESH the finite difference mesh may be either uniform, or user-defined in both
and
directions.
4 References
Hull J (1989) Options, Futures and Other Derivative Securities Prentice–Hall
Wilmott P, Howison S and Dewynne J (1995) The Mathematics of Financial Derivatives Cambridge University Press
5 Parameters
- 1: KOPT – INTEGERInput
On entry: specifies the kind of option to be valued.
- A European call option.
- An American call option.
- A European put option.
- An American put option.
Constraint:
, , or .
- 2: X – REAL (KIND=nag_wp)Input
On entry: the exercise price .
- 3: MESH – CHARACTER(1)Input
On entry: indicates the type of finite difference mesh to be used:
- Uniform mesh.
- Custom mesh supplied by you.
Constraint:
or .
- 4: NS – INTEGERInput
On entry: the number of stock prices to be used in the finite difference mesh.
Constraint:
.
- 5: S(NS) – REAL (KIND=nag_wp) arrayInput/Output
On entry: if
,
must contain the
th stock price in the mesh, for
. These values should be in increasing order, with
and
.
If , must be set to and to , but need not be initialized, as they will be set internally by the routine in order to define a uniform mesh.
On exit: if
, the elements of
S define a uniform mesh over
.
If
, the elements of
S are unchanged.
Constraints:
- if , and , for ;
- if , .
- 6: NT – INTEGERInput
On entry: the number of time-steps to be used in the finite difference method.
Constraint:
.
- 7: T(NT) – REAL (KIND=nag_wp) arrayInput/Output
On entry: if
then
must contain the
th time in the mesh, for
. These values should be in increasing order, with
and
.
If then must be set to and to , but need not be initialized, as they will be set internally by the routine in order to define a uniform mesh.
On exit: if
, the elements of
T define a uniform mesh over
.
If
, the elements of
T are unchanged.
Constraints:
- if , and , for ;
- if , .
- 8: TDPAR() – LOGICAL arrayInput
On entry: specifies whether or not various parameters are time-dependent. More precisely, is time-dependent if and constant otherwise. Similarly, specifies whether is time-dependent and specifies whether is time-dependent.
- 9: R() – REAL (KIND=nag_wp) arrayInput
-
Note: the dimension of the array
R
must be at least
if
, and at least
otherwise.
On entry: if
then
must contain the value of the risk-free interest rate
at the
th time in the mesh, for
.
If then must contain the constant value of the risk-free interest rate . The remaining elements need not be set.
- 10: Q() – REAL (KIND=nag_wp) arrayInput
-
Note: the dimension of the array
Q
must be at least
if
, and at least
otherwise.
On entry: if
then
must contain the value of the continuous dividend
at the
th time in the mesh, for
.
If then must contain the constant value of the continuous dividend . The remaining elements need not be set.
- 11: SIGMA() – REAL (KIND=nag_wp) arrayInput
-
Note: the dimension of the array
SIGMA
must be at least
if
, and at least
otherwise.
On entry: if
then
must contain the value of the volatility
at the
th time in the mesh, for
.
If then must contain the constant value of the volatility . The remaining elements need not be set.
- 12: ALPHA – REAL (KIND=nag_wp)Input
On entry: the value of
to be used in the time-stepping scheme. Typical values include:
- Explicit forward Euler scheme.
- Implicit Crank–Nicolson scheme.
- Implicit backward Euler scheme.
The value gives second-order accuracy in time. Values greater than give unconditional stability. Since is at the limit of unconditional stability this value does not damp oscillations.
Suggested value:
.
Constraint:
.
- 13: NTKEEP – INTEGERInput
On entry: the number of solutions to be stored in the time direction. The routine calculates the solution backwards from
to
at all times in the mesh. These time solutions and the corresponding Greeks will be stored at times
, for
, in the arrays
F,
THETA,
DELTA,
GAMMA,
LAMBDA and
RHO. Other time solutions will be discarded. To store all time solutions set
.
Constraint:
.
- 14: F(LDF,NTKEEP) – REAL (KIND=nag_wp) arrayOutput
On exit: , for and , contains the value of the option at the th mesh point at time .
- 15: THETA(LDF,NTKEEP) – REAL (KIND=nag_wp) arrayOutput
- 16: DELTA(LDF,NTKEEP) – REAL (KIND=nag_wp) arrayOutput
- 17: GAMMA(LDF,NTKEEP) – REAL (KIND=nag_wp) arrayOutput
- 18: LAMBDA(LDF,NTKEEP) – REAL (KIND=nag_wp) arrayOutput
- 19: RHO(LDF,NTKEEP) – REAL (KIND=nag_wp) arrayOutput
On exit: the values of various Greeks at the
th mesh point
at time
, as follows:
- 20: LDF – INTEGERInput
On entry: the first dimension of the arrays
F,
THETA,
DELTA,
GAMMA,
LAMBDA and
RHO as declared in the (sub)program from which D03NCF is called.
Constraint:
.
- 21: WORK() – REAL (KIND=nag_wp) arrayWorkspace
- 22: IWORK(NS) – INTEGER arrayWorkspace
- 23: IFAIL – INTEGERInput/Output
-
On entry:
IFAIL must be set to
,
. If you are unfamiliar with this parameter you should refer to
Section 3.3 in the Essential Introduction for details.
For environments where it might be inappropriate to halt program execution when an error is detected, the value
is recommended. If the output of error messages is undesirable, then the value
is recommended. Otherwise, if you are not familiar with this parameter, the recommended value is
.
When the value is used it is essential to test the value of IFAIL on exit.
On exit:
unless the routine detects an error or a warning has been flagged (see
Section 6).
6 Error Indicators and Warnings
If on entry
or
, explanatory error messages are output on the current error message unit (as defined by
X04AAF).
Errors or warnings detected by the routine:
On entry, | , |
or | , |
or | or , |
or | , |
or | , |
or | , |
or | , |
or | , |
or | , |
or | , |
or | , |
or | . |
and the constraints:
are violated. Thus the end points of the uniform mesh are not in order.
and the constraints:
- , for ,
- , for
are violated. Thus the mesh points are not in order.
7 Accuracy
The accuracy of the solution
and the various derivatives returned by the routine is dependent on the values of
NS and
NT supplied, the distribution of the mesh points, and the value of
ALPHA chosen. For most choices of
ALPHA the solution has a truncation error which is second-order accurate in
and first order accurate in
. For
the truncation error is also second-order accurate in
.
The simplest approach to improving the accuracy is to increase the values of both
NS and
NT.
8.1 Timing
Each time-step requires the construction and solution of a tridiagonal system of linear equations. To calculate each of the derivatives
LAMBDA and
RHO requires a repetition of the entire solution process. The time taken for a call to the routine is therefore proportional to
.
8.2 Algorithmic Details
D03NCF solves equation
(1) using a finite difference method. The solution is computed backwards in time from
to
using a
scheme, which is implicit for all nonzero values of
, and is unconditionally stable for values of
. For each time-step a tridiagonal system is constructed and solved to obtain the solution at the earlier time. For the explicit scheme (
) this tridiagonal system degenerates to a diagonal matrix and is solved trivially. For American options the solution at each time-step is inspected to check whether early exercise is beneficial, and amended accordingly.
To compute the arrays
LAMBDA and
RHO, which are derivatives of the stock value
with respect to the problem parameters
and
respectively, the entire solution process is repeated with perturbed values of these parameters.
9 Example
This example, taken from
Hull (1989), solves the one-dimensional Black–Scholes equation for valuation of a
-month American put option on a non-dividend-paying stock with an exercise price of $50. The risk-free interest rate is 10% per annum, and the stock volatility is 40% per annum.
A fully implicit backward Euler scheme is used, with a mesh of stock price intervals and time intervals.
9.1 Program Text
Program Text (d03ncfe.f90)
9.2 Program Data
Program Data (d03ncfe.d)
9.3 Program Results
Program Results (d03ncfe.r)