NAG FL Interface
d03ncf (dim1_​blackscholes_​fd)

1 Purpose

d03ncf solves the Black–Scholes equation for financial option pricing using a finite difference scheme.

2 Specification

Fortran Interface
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, Intent (In) :: kopt, ns, nt, ntkeep, ldf
Integer, Intent (Inout) :: ifail
Integer, Intent (Out) :: iwork(ns)
Real (Kind=nag_wp), Intent (In) :: x, r(*), q(*), sigma(*), alpha
Real (Kind=nag_wp), Intent (Inout) :: s(ns), t(nt), f(ldf,ntkeep), theta(ldf,ntkeep), delta(ldf,ntkeep), gamma(ldf,ntkeep), lambda(ldf,ntkeep), rho(ldf,ntkeep)
Real (Kind=nag_wp), Intent (Out) :: work(4*ns)
Logical, Intent (In) :: tdpar(3)
Character (1), Intent (In) :: mesh
C Header Interface
#include <nag.h>
void  d03ncf_ (const Integer *kopt, const double *x, const char *mesh, const Integer *ns, double s[], const Integer *nt, double t[], const logical tdpar[], const double r[], const double q[], const double sigma[], const double *alpha, const Integer *ntkeep, double f[], double theta[], double delta[], double gamma[], double lambda[], double rho[], const Integer *ldf, double work[], Integer iwork[], Integer *ifail, const Charlen length_mesh)
The routine may be called by the names d03ncf or nagf_pde_dim1_blackscholes_fd.

3 Description

d03ncf solves the Black–Scholes equation (see Hull (1989) and Wilmott et al. (1995))
f t +r-qS f S +σ2S22 2f S2 =rf (1)
Smin<S<Smax,  tmin<t<tmax, (2)
for the value f of a European or American, put or call stock option, with exercise price X. In equation (1) t is time, S is the stock price, r is the risk free interest rate, q is the continuous dividend, and σ is the stock volatility. According to the values in the array tdpar, the arguments r, q 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 alpha=0.0 and implicit for nonzero values of alpha. Second order time accuracy can be obtained by setting alpha=0.5. According to the value of the argument mesh the finite difference mesh may be either uniform, or user-defined in both S and t 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 Arguments

1: kopt Integer Input
On entry: specifies the kind of option to be valued.
kopt=1
A European call option.
kopt=2
An American call option.
kopt=3
A European put option.
kopt=4
An American put option.
Constraint: kopt=1, 2, 3 or 4.
2: x Real (Kind=nag_wp) Input
On entry: the exercise price X.
3: mesh Character(1) Input
On entry: indicates the type of finite difference mesh to be used:
mesh='U'
Uniform mesh.
mesh='C'
Custom mesh supplied by you.
Constraint: mesh='U' or 'C'.
4: ns Integer Input
On entry: the number of stock prices to be used in the finite difference mesh.
Constraint: ns2.
5: sns Real (Kind=nag_wp) array Input/Output
On entry: if mesh='C', si must contain the ith stock price in the mesh, for i=1,2,,ns. These values should be in increasing order, with s1=Smin and sns=Smax.
If mesh='U', s1 must be set to Smin and sns to Smax, but s2,s3,,sns-1 need not be initialized, as they will be set internally by the routine in order to define a uniform mesh.
On exit: if mesh='U', the elements of s define a uniform mesh over Smin,Smax.
If mesh='C', the elements of s are unchanged.
Constraints:
  • if mesh='C', s10.0 and si<si+1, for i=1,2,,ns-1;
  • if mesh='U', 0.0s1<sns.
6: nt Integer Input
On entry: the number of time-steps to be used in the finite difference method.
Constraint: nt2.
7: tnt Real (Kind=nag_wp) array Input/Output
On entry: if mesh='C' then tj must contain the jth time in the mesh, for j=1,2,,nt. These values should be in increasing order, with t1=tmin and tnt=tmax.
If mesh='U' then t1 must be set to tmin and tnt to tmax, but t2,t3,,tnt-1 need not be initialized, as they will be set internally by the routine in order to define a uniform mesh.
On exit: if mesh='U', the elements of t define a uniform mesh over tmin,tmax.
If mesh='C', the elements of t are unchanged.
Constraints:
  • if mesh='C', t10.0 and tj<tj+1, for j=1,2,,nt-1;
  • if mesh='U', 0.0t1<tnt.
8: tdpar3 Logical array Input
On entry: specifies whether or not various arguments are time-dependent. More precisely, r is time-dependent if tdpar1=.TRUE. and constant otherwise. Similarly, tdpar2 specifies whether q is time-dependent and tdpar3 specifies whether σ is time-dependent.
9: r* Real (Kind=nag_wp) array Input
Note: the dimension of the array r must be at least nt if tdpar1=.TRUE., and at least 1 otherwise.
On entry: if tdpar1=.TRUE. then rj must contain the value of the risk-free interest rate rt at the jth time in the mesh, for j=1,2,,nt.
If tdpar1=.FALSE. then r1 must contain the constant value of the risk-free interest rate r. The remaining elements need not be set.
10: q* Real (Kind=nag_wp) array Input
Note: the dimension of the array q must be at least nt if tdpar2=.TRUE., and at least 1 otherwise.
On entry: if tdpar2=.TRUE. then qj must contain the value of the continuous dividend qt at the jth time in the mesh, for j=1,2,,nt.
If tdpar2=.FALSE. then q1 must contain the constant value of the continuous dividend q. The remaining elements need not be set.
11: sigma* Real (Kind=nag_wp) array Input
Note: the dimension of the array sigma must be at least nt if tdpar3=.TRUE., and at least 1 otherwise.
On entry: if tdpar3=.TRUE. then sigmaj must contain the value of the volatility σt at the jth time in the mesh, for j=1,2,,nt.
If tdpar3=.FALSE. then sigma1 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:
alpha=0.0
Explicit forward Euler scheme.
alpha=0.5
Implicit Crank–Nicolson scheme.
alpha=1.0
Implicit backward Euler scheme.
The value 0.5 gives second-order accuracy in time. Values greater than 0.5 give unconditional stability. Since 0.5 is at the limit of unconditional stability this value does not damp oscillations.
Suggested value: alpha=0.55.
Constraint: 0.0alpha1.0.
13: ntkeep Integer Input
On entry: the number of solutions to be stored in the time direction. The routine calculates the solution backwards from tnt to t1 at all times in the mesh. These time solutions and the corresponding Greeks will be stored at times ti, for i=1,2,,ntkeep, in the arrays f, theta, delta, gamma, lambda and rho. Other time solutions will be discarded. To store all time solutions set ntkeep=nt.
Constraint: 1ntkeepnt.
14: fldfntkeep Real (Kind=nag_wp) array Output
On exit: fij, for i=1,2,,ns and j=1,2,,ntkeep, contains the value f of the option at the ith mesh point si at time tj.
15: thetaldfntkeep Real (Kind=nag_wp) array Output
16: deltaldfntkeep Real (Kind=nag_wp) array Output
17: gammaldfntkeep Real (Kind=nag_wp) array Output
18: lambdaldfntkeep Real (Kind=nag_wp) array Output
19: rholdfntkeep Real (Kind=nag_wp) array Output
On exit: the values of various Greeks at the ith mesh point si at time tj, as follows:
thetaij= f t , deltaij= f S , gammaij= 2f S2 , lambdaij= f σ , rhoij= f r .  
20: ldf Integer Input
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: ldfns.
21: work4×ns Real (Kind=nag_wp) array Workspace
22: iworkns Integer array Workspace
23: ifail Integer Input/Output
On entry: ifail must be set to 0, -1 or 1 to set behaviour on detection of an error; these values have no effect when no error is detected.
A value of 0 causes the printing of an error message and program execution will be halted; otherwise program execution continues. A value of -1 means that an error message is printed while a value of 1 means that it is not.
If halting is not appropriate, the value -1 or 1 is recommended. If message printing is undesirable, then the value 1 is recommended. Otherwise, the value 0 is recommended. When the value -1 or 1 is used it is essential to test the value of ifail on exit.
On exit: ifail=0 unless the routine detects an error or a warning has been flagged (see Section 6).

6 Error Indicators and Warnings

If on entry ifail=0 or -1, explanatory error messages are output on the current error message unit (as defined by x04aaf).
Errors or warnings detected by the routine:
ifail=1
On entry, alpha=value.
Constraint: alpha1.0.
On entry, alpha=value.
Constraint: alpha0.0.
On entry, kopt=value.
Constraint: kopt=1, 2, 3 or 4.
On entry, ldf=value and ns=value.
Constraint: ldfns.
On entry, mesh=value.
Constraint: mesh='U' or 'C'.
On entry, ns=value.
Constraint: ns2.
On entry, nt=value.
Constraint: nt2.
On entry, ntkeep=value.
Constraint: ntkeep1.
On entry, ntkeep=value and nt=value.
Constraint: ntkeepnt.
On entry, s1=value.
Constraint: s10.0.
On entry, t1=value.
Constraint: t10.0.
ifail=2
On entry, sns=value and s1=value.
Constraint: sns>s1.
On entry, tnt=value and t1=value.
Constraint: tnt>t1.
ifail=3
On entry, svalue+1svalue.
Constraint: when mesh='C', si<si+1, for i=1,2,,ns-1.
On entry, tvalue+1tvalue.
Constraint: when mesh='C', ti<ti+1, for i=1,2,,nt-1.
ifail=-99
An unexpected error has been triggered by this routine. Please contact NAG.
See Section 7 in the Introduction to the NAG Library FL Interface for further information.
ifail=-399
Your licence key may have expired or may not have been installed correctly.
See Section 8 in the Introduction to the NAG Library FL Interface for further information.
ifail=-999
Dynamic memory allocation failed.
See Section 9 in the Introduction to the NAG Library FL Interface for further information.

7 Accuracy

The accuracy of the solution f 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 S and first order accurate in t. For alpha=0.5 the truncation error is also second-order accurate in t.
The simplest approach to improving the accuracy is to increase the values of both ns and nt.

8 Parallelism and Performance

d03ncf is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
d03ncf makes calls to BLAS and/or LAPACK routines, which may be threaded within the vendor library used by this implementation. Consult the documentation for the vendor library for further information.
Please consult the X06 Chapter Introduction for information on how to control and interrogate the OpenMP environment used within this routine. Please also consult the Users' Note for your implementation for any additional implementation-specific information.

9 Further Comments

9.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 ns×nt.

9.2 Algorithmic Details

d03ncf solves equation (1) using a finite difference method. The solution is computed backwards in time from tmax to tmin using a λ scheme, which is implicit for all nonzero values of λ, and is unconditionally stable for values of λ>0.5. For each time-step a tridiagonal system is constructed and solved to obtain the solution at the earlier time. For the explicit scheme (λ=0) 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 f with respect to the problem arguments σ and r respectively, the entire solution process is repeated with perturbed values of these arguments.

10 Example

This example, taken from Hull (1989), solves the one-dimensional Black–Scholes equation for valuation of a 5-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 20 stock price intervals and 10 time intervals.

10.1 Program Text

Program Text (d03ncfe.f90)

10.2 Program Data

Program Data (d03ncfe.d)

10.3 Program Results

Program Results (d03ncfe.r)
GnuplotProduced by GNUPLOT 4.6 patchlevel 3 −40 −20 0 20 40 0 20 40 60 80 100 −10 −5 0 5 10 15 Option Values Derivatives Stock Price Example Program Option Values and Derivatives at 5 Months to Maturity option values θ δ γ λ ρ gnuplot_plot_1 gnuplot_plot_2 gnuplot_plot_3 gnuplot_plot_4 gnuplot_plot_5 gnuplot_plot_6
GnuplotProduced by GNUPLOT 4.6 patchlevel 3 −40 −20 0 20 40 0 20 40 60 80 100 −8 −6 −4 −2 0 2 4 6 8 10 12 Option Values Derivatives Stock Price Option Values and Derivatives at 3.5 Months to Maturity option values θ δ γ λ ρ gnuplot_plot_1 gnuplot_plot_2 gnuplot_plot_3 gnuplot_plot_4 gnuplot_plot_5 gnuplot_plot_6