System and method for financial instrument modeling and using Monte Carlo simulation6772136Abstract The software synthesis method and system of the present invention provides a problem solving environment for Monte Carlo simulations (or other concise mathematical description), common in engineering, finance, and science, which automatically transforms a problem description into executable software code. The method and system uses a specification language to support a user's natural description of the geometry and mathematics of the problem and solution strategies. The natural description is concisely expressed using general coordinates and dimensionless parameters, using domain specific keywords as appropriate. The user's problem description is compared with the system's knowledge base to refine the problem--i.e., identifying constraints, applying heuristics and defaults, and applying elaboration rules. The software synthesis method and system uses planning process, computer algebra, and templates to analyze and optimize the problem description, choose and customize data structures, and generate pseudo-code. The pseudo-code is translated into the desired target language source code. The software synthesis system and method therefore provides the ability to describe a problem and possible solution strategies at a high level, and outputs target language code that implements a solution. The software synthesis system and method is particularly useful modeling options where a Monte Carlo simulation is used. Claims What is claimed is: Description BACKGROUND OF THE INVENTION
TABLE 1.4.I.A
Template Declaration for a Conjugate Gradient Solver
Template [ConjugateGradient, SolverTemplate,
SubroutineName[CGSolver],
LocalVars[p[SameAs[y], "search direction vector"],
q[SameAs[y], "image of search direction vector"],
zz[SameAs[y], "conditioned residual"],
alpha[Real, "multiple of search direction vector"],
beta[Real, "orthogonalization constant"],
rho[Real, "temporary variable"],
old[Real, "old value of rho"],
iter[Integer, "iteration counter"],
Templates [Preconditioner],
Code[r = y - svar.xx;
old = r.r;
p = 0;
Until[iter, maxit, StopCode,
zz = r;
commentAbout["Apply preconditioner",
Preconditioner[y -> r. xx -> zz1];
rho = r.zz;
beta = rho / old;
p = zz + beta p;
q = svar.p;
alpha = p.q;
alpha = rho / alpha;
xx = xx + alpha p;
r = r - alpha q;
old = rho]]];
Data declarations in templates are quite different than those of Fortran or C. The System's type annotations, rather than providing fixed storage declarations, allow declarations in terms of configurations of other variables, including those from problem's equations. In this way they resemble the automatic data objects of Fortran 90. Declarations can state that a new array has the same shape as another array, or the same shape but with some dimensions deleted or added, or with dimensionality and range based on the ranges and shapes of other problem or template variables. Templates can be composed just like subroutines in a programming language. The need for a subordinate or "called" template is declared by introducing an identifier, much as local variables are declares. In the pseudo-code, the "call" is represented as a reference to that identifier, signifying only that some class of routine, such as a solver or a preconditioner, is needed at that point. The specific choice of solver or preconditioner is not fixed in the calling template declaration, but rather is a design choice. During synthesis, the template identifier takes on the "value" of an instance of the subordinate template class selected by the specification 'Or system's decision mechanism. For example, in Table 1.4.1A, the identifier "Preconditioner" is declared. Its appearance in the code indicates where to insert the preconditioner call. By default, the System will opt for a no-action preconditioner for the conjugate gradient solver, but users may override this choice by naming any suitably configured preconditioner template. The System's knowledge base represents each template as a class in a hierarchy of template classes. This hierarchy allows for a separation of concern in template declarations as well as defining families of options for various -kinds of algorithms. Among the major classes of templates are solvers, time evolution templates, stepping algorithms, and inverse problem or test wrappers (e.g., for automatically generating convergence rate testing code). Using normal object inheritance, aspects like variables or template place holders defined high in the hierarchy are known in all the subclasses and need not be re-declared. For example, in Table 1.4.1 A the StopCode declaration is inherited by the ConjugateGradient solver from its super-class, Solver Template.
TABLE 1.4.1.B
Class Hierarchy for Evolution Templates
EvolutionTemplate
SteadyState
TimeDependent
TimeDependentEvolve
Movie
Evolve
Motion
DiscreteEventsGeneral
TimeDependentFixed
TestingTemplate
ConvergenceRateTest
EvolveFixed
MotionFixed
Names for both variables and template place-holders can also be drawn from other templates. Mechanisms are available to refer to the variables that exist in the template's "caller" or in the System's problem state, and to do a form of parameter passing from templates to their subordinate templates. A template declaration can also provide an indexing scheme, i.e., a function that transforms the generic indexing operations in the code declaration to more concrete operations appropriate for the synthesized data structures. A template declaration can include certain programming directives. For example, it may specify whether the System should expand the template in line or encapsulate it as a named subroutine.
TABLE 1.4.1.C
Pseudo-Code for PreConditioned ConjugateGradient
seq[assign[r1, g1 - Sal . xx1],
assign[old1, r1 , r1], assign[p1, 0],
seq[assign[iterl, 1],
seq[if[toll < 8*eps1,
print["Warning:Linear tolerance is too small.",
toli, "< 8*", eps1, "."], seq [ ]],
comment ["Compute the static stopping criteria information"],
commentAbout ["The norm of the rignt hand side",
assign[normg1, norm[g1]]],
. . . ]]]
1.4.2 Template Translation The System contains a rich collection of templates written by its developers and submitted to the System via a translator that builds them into the System's object and rule knowledge base. The same translator will be available to users, so that if the templates they need are not already in the System, they may declare them and have them incorporated as first class objects in the synthesis environment. To ensure all the pieces fit together in a sensible manner, the language supports good syntax-checking and type-checking, both at template declaration time and at synthesis time. When a template is declared, the template translator automatically creates an attribute in the template class corresponding to each variable that the template can access and each subordinate template place holder: it also creates methods for filling in the values of these attributes. 1.4.3 Templates During Synthesis During synthesis, the System creates an instance of the selected template class and fills in its attributes as the problem demands. For the highest level components of the algorithm, such as the evolution template and the solver template, the System decides what templates to use as the synthesis process moves from mathematics into algorithm construction. Next, the System fills in the names of the specific variables that can be accessed by extracting them from the problem state or from calling templates, or, in the case of local variables, by defining variables of the appropriate kind and shape. For the attributes corresponding to template place-holders, the methods instantiate the new templates, and the process continues recursively. In the next step of template processing, the System builds pseudo-code for the various template instances. The work is done by methods derived by the template translator from the code part of the template declaration. The first step is fill in the pseudo-code from the declaration with the variables and code for subordinate templates previously derived. Thus, where a subordinate template is called, the template code (see below) for that template is inserted, either in-line or as a subroutine call. The pseudo-code for the conjugate gradient solver in the diffusion example appears in Table 1.4.3. Next, the indexing function, if present, is applied to the pseudo-code. This function transforms the abstract indexing given in the pseudo-code to a concrete form appropriate for the shape of the data structures for the particular variables in the code. No such transformation is specified for the conjugate gradient solver in our example. call[CGSolver, seps1, f, g1, iMax, jMax, SA1] Table 1.4.3: Template Code for Calling ConjugateGradient Finally, if the template is to be cast into a subroutine, the System-encapsulates its code as a subroutine object with appropriate variable declarations and replaces the pseudo-code with a call to the subroutine. This call will be incorporated into the code for the calling templates. If the code is to be in line, no action is needed. The resulting code for calling the conjugate gradient solver in the diffuision example appears in Table 1.4.3. 1.4.4 Data Structure and Operations The ability to generate and customize data structures is one technique that makes the use of flexible abstract algorithm templates feasible. Although the templates declare some minimal constraints on the type of date required, they do not specify the concrete data structures to be used. The System generates these in a manner that customizes them to the specifics of the problem description to minimize the space used by the data structures and the time required to operate on the data structures. The use of mathematical representations of the data structures and their processing by computer algebra both frees the System from having to store libraries of cross products of how to implement all operations on all combinations of data structures and allows an infinite number of representations to be generated based on the specific problem being solved. Although object-oriented assembly allows data structures to be separated from algorithms, component generation goes further in allowing an infinite number of representations to be customized for the specific problem being solved. For example, after equations are transformed into discrete representations, the "stencil"--the number of neighboring points (which can differ in each direction) in an approximation--is used to determine a specialized array representation. The operations are initially represented with sums and products and enumerations over the full array, with the special array represented with conditional constricts. As a much simpler example, a diagonal array f[i,j] is represented by if[i==j, f [i], 0]. Operations are then simplified using computer algebra. For example, the sum of two diagonal arrays, is simplified from doPar[doPar[f[i,j]+g[i,j], <iRange>], <jRange>] via doPar[doPar(if[i==j, f[i], 0]+if[x==y, g'[i], 0], <iRange>], <jRange>] to dopar[f[i]+g'[i], <iRange>] With multi-dimensional stencil arrays and more complex operations such as transpose, the representations and optimizations are much more complicated. A large set of rules are required to simplify expressions involving conditionals, parallel and sequential statements, parallel and sequential loops, comments, and arithmetic and logical operations. This situation-based derivation avoids the necessity for a storing a large library of specialized components to handle the cross product of operations and data structures. 2. Diffusion Example 2.1 Diffusion Problem Specification If a Cartesian coordinate system, with x and y as the spatial variables and t as the time variable is used to describe the problem, then textbooks give the simple diffusion equation as
TABLE 2.1
Problem Specification in Cartesian Coordinates
##EQU1## (Eqn. 1)
(* Geometry *)
Region [0<=x<=2&&0<=y<=3&&0<=t<=1/2, (Eqn. 2)
Cartesian [{x,y}, t]];
(* The PDE *)
When [Interior, (Eqn. 3)
der [f .multidot.t] = = der [f,{x,2}] + der [f, {y,2}];
SolveFor [f]];
(* The BCs*)
When [Boundary, f = = 0]; (Eqn. 4)
(* The IC *)
When [Initial Value, f = = 50 Sin [Pi x/2] Sin [Pi y/3]]; (Eqn. 5)
(* Evolution Algorithm *)
Movie [frames = = 11]; (Eqn. 6)
(* Performance *)
RelativeErrorTolerance [.01];
(* Runtime Interfaces *)
Output [f, Labelled];
Dimensionless;
where the solution f=f(x, y, t) is a function of the space and time variables. If the physical region in x-y-space is chosen where the PDE is to be solved as a simple rectangle with the time interval chosen to be the infinite interval starting at t=0: 0.ltoreq.x.ltoreq.2, 0.ltoreq.y.ltoreq.3, 0.ltoreq.t.ltoreq..infin., (Eqn. 7) The boundary of this region has four pieces. A Dirichlet boundary condition is chosen for the full boundary, that is, the solution is set to zero on each piece of the boundary: x=0, 0.ltoreq.y.ltoreq.3, t.gtoreq.0 f(0, y, t)=0, x=2, 0.ltoreq.y.ltoreq.3, t.gtoreq.0 f(2, y, t)=0, y=0, 0.ltoreq.x.ltoreq.2, t.gtoreq.0 f(x, 0, t)=0, y=3, 0.ltoreq.y.ltoreq.2, t.gtoreq.0 f(x, 3, t)=0, (Eqn. 8) The initial condition, at t=0 and for 0.ltoreq.y.ltoreq.3 and 0.ltoreq.x.ltoreq.2, is chosen to be ##EQU2## It is known that the solution of this problem decays rapidly to zero, so the final value of time t.sub.f to be used in the program can be set rather low, say t.sub.f =1/2. It is desirable to create a movie of the solution; that is, a sequence of three-dimensional plots of the solution f (x, y, t) at equally spaced times. A problem specification that can be used to create a FORTRAN code to solve this problem is given in Table 2.1. This file contains comments delineated as (* comment *), which may be omitted. The statements in the problem specification given in Table 2.1 are direct translations of the mathematical statements into specifications in the problem specification language. The first argument to Region is a Boolean combination, using the logical and operator &&, of inequalities, while the remaining two arguments are optional. The Boolean argument is used to describe the space-time region; the first optional argument specifies the coordinate system and the second gives the name of the region. The region must be specified using less than or equals (<=) or greater than or equals (>=) rather than strict inequalities. Cartesian takes one or two arguments. The first is a list of the spatial coordinates, of which there may be one, two, or three. The second argument is the name of the time variable. This specification produces a Cartesian coordinate system with axes labeled by the members of the first argument. The PDE, that is, the diffusion or field equation, is given using standard mathematical operators from Mathematica as well as der, which replaces the Mathematica D. As with D, der specifies some derivative of its first argument; the remaining arguments are lists of two items where the first is the variable with respect to which the derivative is to be taken and the second is the order of the derivative, for example: ##EQU3## The == is the standard mathematical equality which coincides with Mathematica's notation. The use of der rather than the Mathematica D points out one of the shortcomings of most computer algebra systems for tasks such as discretizing differential equations. Mathematica's D is designed to compute derivatives. Most of the time in this application, the derivatives are formal, that is, derivatives of an abstract function. Typically the function is to be computed numerically, so the properties of the D operator are not useful, and in fact, frequently produce results exactly the opposite of what is needed. For example, it usually not appropriate to apply the chain rule. Also, in coordinate free notation it is desirable to omit the dependencies, which would normally cause D to return 0. When [<region>, <eqn>] specifies an equation and a region over which the equation holds. It takes two arguments: the first is a Boolean that describes pieces of the space-time region, or a name of a previously described region, while the second is a partial differential equation, a boundary condition, or an initial condition or a name of a previously described equation. For more complex problems, systems of equations can replace a single equation. The specifications Initial Value, Interior, and Boundary are shorthand, respectively, for the physical region at the initial time, the interior of the full space-time region, and all of the boundary of the full space-time region. For example, if the region is given the name R as in Region [0<=x<=2 && 0<=y<=3 && 0<=t<=2, Cartesian [0{x, y}, t], R] then Interior means Interior [R], which is equivalent to the Region 0<x<2 && 0<y<3 && 0<t<2. The boundaries and the equations for the boundary conditions can also be given explicitly; in this example, an equivalent specification to When [Boundary, f=0] would be When [x==0, f>=0]; When [x==2, f==0]; When [y==0, f==0]; When [y==3, f==0]; SolveFor specifies what variables to solve for; it can be omitted if the solution variables are obvious. Movie says to integrate an evolution problem with output at the ends of equal time intervals, the number of which are given by frames==11. A specification of the results to output is given by the expressions in Output. In this case, the values for f in the full region, are written into a file that can be used to plot the solution. The argument Labelled specifies that the values of time and the spatial points where f is computed are included in the file. Far more detailed control over the output is available. The RelativeErrorTolerance gives a heuristic tolerance for the relative error. The System does not guarantee that this tolerance will be met. The Dimensionless specification tells the System the problem is dimensionless, thereby suppressing checks for the consistency of physical units in equations. For straightforward problems, this is all that needs to be specified. There are many other defaults that the System will use to translate the problem specification given in Table 2.1 to a numerical code. For example, the numerical code produced is, by default, is FORTRAN-77, but this can be changed to C for example. Also, the System will select what it considers appropriate Discretization methods and solvers if the specification does not indicate what to use. It is easy to modify the problem specification in Table 2.1 to solve the convection diffusion equation. ##EQU4## 2.2 A Coordinate-free Description of Diffusion The equation specification given in Table 2.1 is only useful in Cartesian coordinates, so if a problem were given on a circular wedge, a completely new problem specification would be needed and the diffusion equation would have to be manually transformed into polar coordinates. In fact, another significant use of computer algebra in the System is writing partial differential equations in Mathematica's thirteen standard orthogonal coordinate systems and in general coordinates so that automatic grid-generation techniques can be used. The best approach, therefore, is to write the differential equation and as many of the initial and boundary conditions as possible in a coordinate-free notation and then let the System figure out how to write them in the desired coordinate system. In coordinate-free notation, the PDE is expressed mathematical as ##EQU5## where a parameter or parameter variable .sigma. has been added, and where .gradient..sup.2 is the Laplacian; in two-dimensional Cartesian coordinates, the Laplacian is ##EQU6## The physical region is still a rectangle, but now we give the length of the sides with the parameters .alpha. and .beta. is: 0.ltoreq.x.ltoreq..alpha., 0.ltoreq.y.ltoreq..beta., 0.ltoreq.t.ltoreq..infin. (Eqn. 14) Also, one of the zero Dirichlet boundary conditions is changed to the parametric value A: f(0,y,t)=0, f(.alpha.,y,t)=0, f(x,0,t)=0 f(x,.beta.,t)=A (Eqn. 15) while the initial condition is chosen to be ##EQU7## where .beta. is a parameter in the IC. (In this case, if A.noteq.0, the boundary condition is not consistent with the initial condition. However, this is not a problem for the diffusion equation.) Again, the problem is to create a movie of the solution using the values for the parameters that were used in the previous problem (see Table 2.1), but with A=100. The System problem specification for solving this problem is given in Table 2.2. The Incline specification tells the System to replace the listed variables by their values in the code (all instances of A will have 100 substituted, and so on). Computer algebra is used to propagate the effects of constant substitution; entire equations are eliminated if appropriate. The expression Laplacian is the notation used by the System for the Laplacian. The System also knows about the divergence, gradient, and curl. Note that parameters can be included in the equations as variables whose values can then be set in other equations. When a formula like sigma==0.01 is used, it typically results in code for setting the value of the parameter sigma occurring early in the numerical program, but parameter values can also be read from the terminal and from files. Because .sigma. is small, t.sub.f, (tfinal) should be taken large. Use of parameters facilitates changing problems quickly. Note that the form When [Boundary, f==0] acts the default boundary values for f to zero on the full boundary. This is then overridden on the piece of the boundary with y=.beta. by using the form beginning When [y==beta, . . . ]. The specification OneFile given in the Output specification will cause all of the data to be written to a single file. The user has not specified the numerical algorithm to be used in solving this problem. Therefore, the System will use its defaults, and employ a simple forward difference operator in time (first-order accurate and explicit), and central difference operators in space (second-order accurate). The resulting algorithm is simple, not requiring the solution of a large linear system). However, it has the well known disadvantage for diffusion equations that the time step must be deceased in proportion to the square of the spatial grid size.
TABLE 2.2
Coordinate-Free Problem Specification
(*Geometry*)
Region [0<=x<=alpha && 0<=y<beta && 0<=t<=tfinal,
Cartesian [{x,y},t]]; (Eqn. 17)
alpha == 2.0; beta == 3.0; tfinal == 10.0
(*The PDE*)
When [Interior, der [f,t] == sigma laplacian]f]; (Eqn. 18)
sigma == .1;SolveFor[f]];
(*Default Boundary Condition*)
When [Boundary,f == 0]; (Eqn. 19)
(*Specific Boundary Condition*)
When [y == beta, f == A; A == 100]; (Eqn. 20)
(*Initial Condition*)
When [Initial Value,
f == .beta. Sin [Pi x/2] Sin [Pi y/3];B == 50]; (Eqn. 21)
(*Evolution Algorithm*)
Movie [ frames == 5];
(*Performance*)
RelativeErrorTolerance [.1];
(*Runtime Interfaces*)
Incline [ {alpha, beta, tfinal, sigma,A,B}];
Output [f, Labelled, OneFile];
Dimensionless;
2.3 Variable Coefficients and Heterogeneous Material In addition to specifying numerical values for parameters, parameters can also be functions of other variables. In the previous examples, the independent variables are x, y, and t, while f is the dependent variable which is to be solved for as a function of the independent variables. In the previous example, the variables .sigma.,.beta.,.beta.,t.sub.f, A, and B are parameters. In this example, a heterogeneous material property is defined to illustrate how to define more complex parameters. Heterogeneous materials can be modeled using coefficients that are defined by functions with a jump. For example, the coefficient sigma in the previous problem, could be defined with jump functions. ##EQU8## However, because the coefficient c is not constant, the differential equation must be written in terms of the divergence and gradient. To solve this problem, the specifications following the comment (* The PDE *) in the problem specification given in Table 2.2.1 are replaced by (* The PDE *) When [Interior, der [f,t]==div [sigma grad[f]] sigma==if[x>1 && y>2, 1, 1/100]]; Note that the "if" statement has three arguments: the first is a Boolean; if the first argument is true, then the value of the second argument is returned; while if the first argument is false, then the value of the third argument is returned. Everything else in the problem specification is unchanged from Table 2.2.1. The fact that the coefficient sigma has a jump has a strong impact on the type of algorithm that should be used to solve the problem. The jump-coefficient effect can be achieved in other ways, such as specifying multiple regions, each with their own equations, or specifying the values of coefficients in material property files. 2.4 Specifying the Solution Procedure In the specification for the heat diffusion problem, Table 2.2.1, the following specifications were used to control the solution procedure and consequently the kind of numerical code that the System generates: (* Evolution Algorithm *) Movie [frames==11]; (* Performance *) RelativeErrorTolerance [0.01]; (* Runtime Interfaces *) Output [f, Labelled, OneFile]; These specifications can be modified or augmented by additional information. A main feature of a code for the solution of an initial-boundary value problem is the overall time-stepping loop which is given by an evolution algorithm. The evolution algorithms are recorded in evolution templates. The next issue in designing a solution algorithm is the time discretization and a spatial discretization. The main types of time discretizations are explicit and implicit, and the typical space discretization is central. Implicit discretizations require the use of a linear equation solver. These choices complete the overall algorithm. Other important details can also be specified. The performance criteria control the accuracy and amount of work required to solve the problem. The main point is to choose either a relative error tolerance for the solution or to directly specify the grid size. There are currently two programming languages available: Fortran-77 and C. The run-time interfaces include input and output. 2.4.1 Evolution Algorithms The evolution algorithm controls the structure of the overall time stepping loop. Some of the critical features are how to choose the time step, in particular, whether to use a fixed or variable time step, and where to print the results of the computation. Changing the Movie specification to Evolve will cause output to be generated only at the beginning and end of the time interval. There are several other evolution templates available. 2.4.2 Time Discretizations For diffusion equations, the default discretization is the simple explicit or forward Euler discretization. This is chosen because it is simple to implement and effective for small problems, but this has the disadvantage that a large number of time steps is required. Diffusion equations are more typically discretized using implicit methods. This produces an additional problem that a linear solver is needed for the resulting system of equations. The System will choose one if none is given in the specification. The Crank-Nicholson time discretization is commonly used and this can be specified by adding CrankNicholson: to the specification file. Now there is no stability constraint, but accuracy requires that the number of time steps be proportional to the number of spatial steps. The default linear solver is Gauss-Siedel (SOR with the relaxation parameter .omega.=1). In this simple problem, a tri-diagonal solver is better choice. This choice can forced by adding TriDiagonal; to the specification file. 2.4.3 Spatial Discretizations For diffusion equations, the default spatial discretizations are second-order central differences, but there are many other discretizations available. Here is another point where computer algebra plays a critical role. The discretizations are given as Mathematica replacement rules that apply to a system of equations, or just to one equation, or if necessary, only to a particular term in an equation. Many discretizations are developed by using the series expansion capabilities of computer algebra. Also, after the discretization is complete, the System uses Mathematica functions to simplify the discretized equations and collect the coefficients of the discretized unknown function, which are known as the stencil coefficients for the discretized equation. In addition to applying the built in discretization rules, users can define their own replacement rules for special discretizations. 2.4.4 Performance Criteria Performance criteria are related to the choice of grid sizes and accuracy tolerances and are the major factors in the amount of work required to solve the problem. The System has some simple heuristics to choose a spatial-discretization size given a relative error tolerance as in RelativeErrorTolerance [0.001]; Note that using the specification RelativeErrorToicrance [retol]; without specifying a value for retol, will require the relative error tolerance to be input at the terminal when the numerical code is run. Given the spatial discretization size, the System then determines a time-discretization size. However, the grid sizes can be specified directly, for example, to test several grid sizes, in which case grid size values can be read from the terminal. However, care should be taken not to violate any stability constraints. For example, the relative error tolerance specification can be replaced by nMax==4*iMax**2; GridSize [{{iMax, iMax, nMax}}]; Note that no value is given for iMax, so when the Fortran code is run, it will request a value for iMax, the spatial discretization size. 2.4.5 Programming Languages There are currently two numerical programming languages available: the default Fortran-77 and C, which are specified by Fortran; or Clang; The numerical precision in the generated codes can be controlled by using Single; which is the default, or Double; The System can be extended to other programming languages. 2.4.6 Runtime Interfaces Runtime interfaces refer to how input and output are handled, whether or not subroutines are to be generated as separate code or are to be in-lined, and how the code generated by the System will interface with other codes. The most obvious way to initialize a variable is to provide an equation for it in the specification, as was noted the examples given previously. In most programming languages, variables can be scalar, typically real or complex, or any array of real or complex numbers. The System generates code appropriate for the type of variable. For example, scalar variables that are mentioned in a special dictation but not given a value are, by default, read from the terminal, while array variables are, by default, read from a file. The procedure for scalar variables is illustrated in the examples. However, it is also possible to output the gradient of the solution or the solution and its gradient, for example, the specification Output [{f, grad[f]}, "Data", Labelled]; will cause t, x, y, f .differential.f/.differential.x, and .differential.f/.differential.y to be written to some output files. The string "Data" determines the first part of file names and the data for each value of t is written to a different file because the specification OneFile was omitted. The specification given in FIG. 2.3 can be modified to put the negative of the gradient of the solution using Output [{-grad [f]}, "Data", Labelled]; The gradient can be plotted as in FIG. 4 using the Mathematica PlotField package. Note that the plots in FIG. 4 clearly show that the initial heat distribution 3. Financial Instruments 3.1 Background New types of exotic options and new twists on familiar ones proliferate as customers demand instruments structured to specific requirements. Models for these new instruments must be generated quickly and inexpensively. Although some analytic formulas exist, most are based on assumptions, e.g. constant volatilities or continuous sampling of spot prices, that may lead to unacceptable errors when applied to practical pricing problems. Most often, analysts rely on such numerical methods as binomial trees, Monte Carlo simulation, or finite difference solution of the governing partial differential equations (PDEs). The current method of generating computational pricing models, like scientific modeling software in general, is conventional programming: a tedious cycle of coding, debugging, testing, and re-coding, quite ill-suited to the present demands. The System and method-of the present invention provides a process of automatically transforming very high level specifications that mirror the mathematical description of the problem into efficient code in a conventional programming language such as C or Fortran, thereby eliminating the programming step completely. The System is well suited for problem in modeling financial instruments. The System synthesizes finite difference codes for general PDES, and produces pricing models for exotic options. The System displays numerical results from the generated codes and compares them to exact and approximate analytic results and to Monte Carlo results. A specification for a complex option might occupy less than a page and yet generate (in less than an hour) a commented, modular C code of several thousand lines, yielding several orders of magnitude increase in software productivity. The System can apply finite difference techniques to complex option pricing problems that have been traditionally attacked with Monte Carlo or tree methods. These latter methods are popular largely due to the perception that they are easy to program. When specifications are made at a high level and programming considerations eliminated, the various methods can be judged on speed and accuracy alone. In many cases of practical interest, finite difference methods have significant performance advantages over other methods. 3.2 Pricing Options with Finite Difference Methods The use of finite difference methods for solving partial differential equations pre-dates the modem computer. A huge body of experience exists in such diverse applications as computational fluid dynamics, electromagnetics, and seismic wave propagation. Convergence rates are generally rapid and well understood. The Black-Scholes equation commonly used for options valuation and its generalizations are ideal candidates for finite difference solution since they are generally linear and because they generally contain dominant diffusive terms that lead to smooth solutions. (Although linearity simplifies the finite difference solution of PDEs, it is certainly not a requirement.) Tree methods have an intuitive appeal, and it is easy to show that trinomial trees are closely related to explicit finite difference methods. However, trees for options involving multiple assets or interest rate securities with arbitrary interest rate process models can be difficult to construct, while the corresponding finite difference models are generally straightforward. It is easy to generate finite difference schemes that are second order accurate, i.e., whose error decreases as the square of the number of time steps. Tree pricing methods are generally first order accurate. With respect to Monte Carlo techniques, finite difference methods have three important advantages, the accurate calculation of hedging parameters, ease of implementing early exercise features, and speed of computation. Because the Black-Scholes equation is so well suited to finite difference solution, the System requires only fairly simple algorithms for most option pricing problems (relative to fluid dynamics, for example). An introduction to finite difference methods for solution of general initial-boundary value problems may be found in (A), while a thorough discussion of finite difference methods specifically for option pricing is found in (B) (all references cited herein incorporated by reference for background). 3.3 Example 1: An American put 3.3.1 Abstract Specification When producing modeling codes, scientists and analysts begin with an abstract specification of the problem to be solved and the techniques to be applied. This specification can often be described to a colleague in a short conversation, yet it conveys virtually all significant aspects of the desired code. Suppose the problem is to generate a finite difference model of American put options. An abstract specification might be: "Solve the Black-Scholes equation in a region with one Cartesian space dimension S, the underlying spot price, and one time dimension t. At the minimum price boundary the option value is the strike price K, and at the maximum PCE boundary the option value is zero. The initial condition is the standard payoff. Constrain the option value to be greater than or equal to the payoff. Use Crank-Nicholson differencing for stability, and a SOR solver to invert the resulting matrix. Write the code in C using double precision arithmetic." The mathematical expression of this abstract specification consists of the Black-Scholes equation (using backward time), pertinent initial and boundary conditions, and constraints: ##EQU9## Despite this mathematical rigor, many details, principally algorithm and programming considerations, have been omitted, to be supplied by common modeling practice. While human translation of such specifications into code are routine, the process is error-fraught and time consuming. 3.3.2 The Program Synthesis Specification
TABLE 3.3.2
Example Specification for American put
Region[0<=S<=SMax && 0<=t<=T, Cartesian[ {S },t]];
When[Interior, CrankNicholson;
der[V,t]== {fraction (1/2)} sigima 2 S 2 der[V, {S,2} ] +
(r-D0) S der[V,S] - r V]];
payoff==Max[K-S,0];
When[min[S], V==K];
When[max[S], V==0];
When[min[t], V==payoff];
Constraint[V >= payoff];
SOR;
TargetLanguage[C];
Double;
The System automatically generates the required numerical code from a specification, shown in Table 3.3.2, that is remarkably similar to the abstract description. The Region (solution domain) is defined in terms of spot price S and backward time t. For an option involving multiple underlying assets, the domain would be multi-dimensional. The System can handle an arbitrary number of spatial dimensions. Users can also specify transformations to generalized coordinates which concentrate grid points near 'specific points such as strike or barrier prices. The Black-Scholes equation is specified in the Interior, i.e., the "not-boundary" part of the region. For entering PDE's in text form, the System uses the notation der[V,S] for .differential..nu./.differential.s. S is the spot price, sigma the volatility, "r" the risk-free interest rate, and D.sub.0 the continuous dividend rate. The boundary condition at the minimum price boundary, S=0, is the usual one, while at the maximum S boundary the option value is simply set to zero, with the understanding that users will set SMax, the maximum value of S, large enough (typically about four or more times the strike price) that the induced error will be negligible. The initial condition for our solution is the usual payoff at expiry for a put, Max[K-S,0l, where K is the strike price, and integrated in backward time from expiry to the present value. 3.3.3 Discretization Specifying a discretization method, CrankNicholson, in the Interior, where it will determine the discretization of the Black-Scholes equation. This implicit differencing method is unconditionally stable and second order accurate in the time step .DELTA.t. Since not specified otherwise, the System will use its default second-order central space differencing in the S direction. Discretization methods for systems of PDEs, individual PDEs, or even selected terms of a PDE may be chosen by the analyst (or by the System if omitted from the specification) from an extensive suite of pre-defined numerical methods, including splitting, Alternating Direction Implicit, and Crank-Nicholson. Alternatively, users may design custom discretizations. For example, in a multidimensional model, such as an outside barrier option, one may want fourth-order differencing in the spot price dimensions and Crank-Nicholson differencing in time, except for mixed derivative terms (those involving the price correlations of the underlying) which are to be second-order in spot price and explicit in time. Such a novel discretization is specifiable in a few lines. 3.3.4 Solvers and Constraints As in this example, most practical numerical algorithms for Black-Scholes and related equations involve implicit differencing for numerical stability, resulting in a linear system of equations to solve. Linear solvers may be selected from about a dozen predefined algorithms known to the System or called from a static library of the user's choosing. The predefined solvers are represented as pseudo-code templates that are expanded incline as the System builds code, taking advantage of knowledge of the specified data representation and properties of the associated matrix. Various preconditioners are also provided. The pseudo-code representations of algorithms in (including solvers) are abstract, devoid of references to specific dimensionality, data structures, equations, or target language, but with notions of parallelism (incorporated for future generation of codes in parallel and distributed languages.) The American feature requires the application of a constraint: the option value cannot be smaller than the payoff. The appropriate formula is simply given as an argument to the Constraint statement. Constraints may be applied inside the selected solver as it iterates toward a solution, or outside the solver after convergence has been achieved. The former approach is better. When the constraint is applied inside a successive over-relaxation solver (the so-called projected SOR method) the result can be shown to converge to the correct solution. In this example, a SOR solver is specified. The constraints are applied inside the solver, and in fact, anywhere in the code that the values of constrained variables are reassigned. The choice of over-relaxation factor, omega, has a critical influence on the performance of SOR solvers. Users may specify solver parameters, e.g. SOR[omega==1.25, maxit==100, tol==tol=10 (-5)]. The optimum value of omega, however, is problem specific and difficult to estimate. The System provides an iterative refinement algorithm for 0,SOR[RefineOmega], which may result in a substantial reduction in the R number of iterations. An optimized value for w is generally more effective in problems that do not have constraints. Finally, the System specifies the precision and target language of the code to be generated, namely Double and C. Alternative choices for target language are Fortran and Fortran90. The ten-line specification in generates (in less than ten minutes on a high-end PC) about three hundred lines commented, modular ANSI standard C code for finite difference valuation of American puts. 3.3.5 Improved Specification
TABLE 3.3.5
Specification for American put, with greeks,
improved initial condition.hz,1/32
Region[0<=S<=SMax && 0<=t<=T, Cartesian[ {S},t]];
When[Interior, CrankNicholson;
BlackScholes1D[V->P, Keywords->{Vega,Rho}]];
payoff == Max[K-S,0];
Constraint[P >= payoff, if[P==payoff, Vega==0;
Rho==0]];
When[min[S], P==K;
Rho==0;
Vega--==0];
When[max[S], P==0;
Rho==0;
Vega==0];
When[min[t], Vega--==0;
Rho==0;
P == if Abs[S-K] > delta[S]/2,
(K-S+delta[S]/2) 2/(2 delta[S]), payoff]];
SOR;
Default[ReadFile["pinit.dat"]];
Greekout[V->P];
TargetLanguage[C];
Double;
Next, the System incorporates some additional features in this model to improve its accuracy and utility. Table 3.3.5 shows an improved specification. The first new feature is the "equation generator" BlackScholes1D, used to specify the Black-Scholes equation. It is convenient to place the definitions of frequently used equations in a separate specification file. Simply mentioning the name of the file expands it as a macro. Typically, there are several forms of the equations of interest, or as in this case, additional equations, which are selectable via Keyword arguments in the macro expansion. Macros may be defined for any sequence of frequently used specification commands, and they may even be arbitrarily nested. Here, the System used a substitution rule, V.fwdarw.P, in expanding the macro, to rename the option value variable from the default V to P. 3.3.6 Hedging Parameters, Equation Generators In order for the American put model to be useful for hedging as well as pricing it must calculate the hedging parameters, i.e. the greeks. The System can use symbolic algebra (e.g. Mathematica) to derive new equations from existing ones. For example, it can differentiate the Black-Scholes equation with respect to any of its parameters to form equations for the greeks. These new equations are then automatically discretized and solved along with the Black-Scholes equation. Boundary and initial conditions for the new quantities must be provided, of course. In Table 3.3.2 the System invokes BlackScholes1D with the Keywords Vega and Rho, directing that two additional equations for these quantities be produced and solved. The greeks, Delta, Gamma, and Theta do not require PDE's that are integrated in time. While the System could generate such equations, it is wasteful to do so. These parameters can be derived from the instantaneous option value by differentiation with respect to S or t. This simple differentiation can be done in the Output commands discussed below.
TABLE 3.3.6
Equation Generate Blackscholes1D
PhysicalMeaning[S,Dollars);
Create the Dependent Variables
Variable[SMin, Scalar, Dollars, "minimum price"];
Variable[SMax, Scalar, Dollars, "maximum price"];
Variable[K, Scalar, DoIlars, "strike price"];
Variable[sigma, Scalar, Time (-1/2), "volatility"];
Variable[r, Scalar, Time (-1), "interest rate"];
Variable(D0, Scalar, Time (-1), "continuous dividend yield"];
Variable[V, Scalar, Dollars, "option value"];
(* Create and name Black-Scholes Equation for option value V *)
Equation[-der[V,t] + (1/2)sigma 2 S 2 der[V, {S,2} ] +
(r-D.sub.0) S der[V,S] - r V==0, ValEqn];
(* Create and name the equation for Vega=der[V,sigma] *)
IsAKeyWord[Vega,
Variable[Vega, Scalar, Dollars Time (1/2), "Vega: der[V,sigma]"];
Equation[Greek[ValEqn, {Vega,der[V,sigma]}], VegaEqn]
];
(* Create and name the equation for Rho=der[V,r] *)
IsAKeyWord[Rho,
Variable[Rho, Scalar, Dollars Time, "Rho: der[V,r]"];
Equation[Greek[ValEqn, {Rho,der[V,r]}], RhoEqn]
];
An equation generator is simply a macro file that encapsulates the definition of an equation, its optional variants, and other related declarations. Table 3.3.6 shows the equation generator file BlackScholes1D.eqn. the System expects coordinate variables to have dimensions of length. The first line of this file simply informs the System that the coordinate variable S has units of Dollars. Next, we give the system some information about the important parameters such as sigma and r, namely their tensor order, their units, e.g., {Dollars}, and a definition string. The tensor order and units are used by the system to check the equations and boundary conditions for tensor order and units consistency. If tensor order and units information are not given, the system uses the defaults Scalar and Dimensionless. The definition strings are incorporated into comments in the generated code. Next, the Black-Scholes equation is specified and given a name, ValEqn, for later reference. If Vega is present in the Keyword list, the variable Vega is created with appropriate units, and a new equation for it, VegaEqn, is generated automatically by operating on ValEqn with the differentiating function Greek. It needn't be specified by hand. Finally, an equation for Rho is created if the KeywordRho is specified. Generalization of this technique to define equations for other hedging parameters is straightforward. Thus, invoking BlackScholes1D with both keywords defines these additional equations: ##EQU10## where .rho. and {character pullout} are Rho and Vega respectively. Returning now to the specification in Table 3.3.6, the System supplies constraints for the hedging parameters. The addition to the Constraint statement simply says that if the option value P has reached its constraint, then it can no longer depend on volatility or interest rate, so Rho=Vega=0. 3.3.7 Improved Initial Condition A small change to the initial condition for P in Table 3.3.5 increases the accuracy and consistency of the model. In the continuous world of the Black-Scholes equation, the initial condition for the option value is just the payoff function we have defined. However, in the discrete world of the finite difference solution, using P[S,0]==payoff generates small errors near the strike price K where payoff has a discontinuity in slope. The final option value depends slightly on whether K falls on a grid point or between grid points. This dependence can be largely eliminated by slightly modifying the initial option value near the strike so that the area under the payoff curve (in the sense of "trapezoidal rule" integration) is conserved delta[S] is the System's text notation for the discretization increment in the S dimension, .DELTA.S. The improved initial condition specification for P uses the conditional statement if[ . . . ] to set P=payoff at all but the single grid point closest to K, where P is modified so that the numerical and analytic areas under the payoff curve are identical. The general form for the if statement is i[test, statements, statements], where the first group of statements generates code for the case that test is true, while the second group of statements generates code for the case that test is False. 3.3.8 Input/Output Finally, the specification of Table 3.3.5 directs the system to handle input and output in specific ways, rather than relying on system defaults. The System sets the default so that any variables such as the strike price K or maximum price SMax, not explicitly defined by an equation will be read from the file pinit.dat. Output can be tailored for convenient entry into most any visualization or post processing package. Because the System reuses the output specifications for the option value and the greeks so often, it can put them in a macro file. The file Greekout in thisexample is such a macro. It contains statements such as .backslash.Output[der[V,S], "delta.out""]+. Placing expressions in Output statements is how it defines and outputs the hedging parameters. FIGS. 5(a) and (b) show the Results for the American put where (a) depicts expiration (dotted curve) and present (solid curve) values, and (b) shows expiration (dotted curve) and present (solid curve) values of, .DELTA. Parameters used: T=1.0, D0=0.025, sigma=0.5, r=0.05, K=1.0, iMax=100 (number of Sgrid points), nMax=20 (number of t grid points), SMax=4.0. FIGS. 6(a) and (b) show the Results for the American put, where Mean (a) and standard deviation (b) of the errors (x 10 4) in the present option value with simple initial condition (dotted curve) and improved initial (solid curve) condition, as the strike K is varied around a grid point. 3.3.9 Numerical Results FIGS. 5(a) and (b) show the initial and final values of P and Delta calculated by the synthesized finite difference code. The American feature has the expected effect, increasing the value of the put. FIGS. 6(a)(b) show the displays the mean error in the final option value for eight slightly different values of the strike K in the range 1-delta[S]/2 K 1+delta[S]/2. The dashed curve is the simple initial condition specified in Table 3.3.1 while the solid curve is the improved initial condition in Table 3.3.2. The mean values of the error are approximately the same, for either initial condition, a few basis points. However, the standard deviation of the error, shown in FIG. 5(b)), is reduced by more than an order of magnitude for the improved initial condition. The small variation of option value with the position of K relative to the finite difference grid is nearly eliminated, falling to a tiny fraction of a basis point. (A very narrow region in which standard deviation is not significantly changed, marks the location of the American free boundary.) Since an analytic solution does not exist for this option, error has been measured relative to an identical calculation that has ten times the number of Sgrid points. 3.3.10 Switches A final enhancement to the put option specification is to incorporate a switch to allow generation of three distinct codes from one specification, shown in Table 3.3.4. The System declares a Boolean variable american. To generate a code for a European option, the system simply types american.fwdarw.False before beginning synthesis. The min[S]boundary condition is changed appropriately, the Constraint code is not generated, and a direct TriDiagonal solver is used instead of the slower, iterative solver SOR. If american.fwdarw.True is defined before synthesis, the previous model is recovered. Finally, if leaving american undefined at synthesis time, a code is generated that contains branches for both European and American options. American then becomes an input variable, read from the default input file. Its value at execution time determines the type of option priced.
TABLE 3.3.10
Specification for American/European put, with greeks
Boolean[american];
Region[0<=S<=SMax && 0<=t<=T, Cartesian {S},t]];
When[Interior, CrankNicholson;
BlackScholes1D[V->P, Keywords->{Vega Rho}]];
payoff == Max[K-S,0];
if[american, Constraint[P >= payoff; if[P==payoff; Vega==0; Rho==0]]];
When[min[S], Rho==0; Vega==0; if[american, P==K, P=32 K Exp[-r t]]];
When[max[S], P==0; Rho==0; Vega==0];
When[min[t], Vega==0; Rho==0;
P == if[Abs[K-S]>0, (K-S+delta[S]/2) 2/(2 delta[S]), payoff]];
If american, SOR, Tridiagonal];
Default[ReadFile["pcinit.dat"]];
Greekout[V->P];
TargetLanguage[C]; Double;
3.4 Example 2 Discrete Knockout Call in Generalized Coordinates Like the American option evaluated in the first example, barrier options are path dependent, but they provide very different numerical challenges. Both Monte Carlo methods and trinomial trees show poor convergence for barrier problems. Several methods have been developed to address these problems. For trees, one obvious step is to ensure that the tree nodes fall exactly on the barrier price. This is not always easy to do, especially in higher dimensions. With finite differences this is a simple matter. However, this example generates generalized grids that concentrate grid points near the barrier. This is accomplished with a few extra lines of specification, but results in vast improvement in accuracy. This example also introduces discrete sampling of path information, e.g., activation of a knockout barrier only at the end of the trading day. Using program synthesis, complex code segments such as those that model the discrete sampling of a barrier, the payment of a discrete dividend or sampling of the spot price for a lookback option are generated from pseudo-code templates. These templates are invoked via simple keyword commands in the specification. 3.4.1 Generalized Coordinates General coordinate transformations are used in the numerical solution of PDE's to increase the density of grid points in critical regions and to allow coordinate systems to conform to complex boundaries. Adaptive grids are generated by time-dependent general coordinate systems to conform that adapt to the evolving solution. The System currently provides for generalized coordinate transformations specified by formulas. Automatic numerical generation of grids with various smoothness and orthogonality properties is under development. The following modes a European knockout call with a discretely sampled barrier using a coordinate transformation to concentrate grid points near the barrier. The mathematical specification of the problem is ##EQU11## where B is the knockout barrier price, and t; is the set of barrier sampling times. The specification is given in Table 3.4.1.
TABLE 3.4.1
Specification for European knockout call in generalized coordinates.
Region [0<=xi<=1 && 0<=tau<=TMax, GeneralCoordinates[xi,tau]];
OldCoordinates[ {{S},t}];
S== if[xi<=xiB,
Smin + alpha c1 (1 -Exp [-A xi]);
Smax - alpha c2 (1 -Exp[A (xi-1)]);
T=tau;
c1= 1+(B-SMin)/alpha;
c2==1+(SMax-B)/alpha; A==Log[c1 c2];
xiB==Log[c1]/A];
payoff--Max[S-K,0];
When[Interior, CrankNicholson; BlackScholes1D[ ]];
When[min[xi].parallel. max[xi], V==0],
When[min[tau], V==payoff];
TriDiagonal;
DiscreteEvents[ Barrier[direction[S],
functions[{if[S>=B, V==0]}],
ReadFile[tsample, "tbar.dat"],
nsample==nbar]];
Default[ReadFile["barinit.dat"]];
Output[V, "V.out", OneFile, Labelled];
TargetLanguage[C]; Double;
The System specifies GeneralCoordinates [xi,tau] as the coordinate system, where xi .xi. and tau (.tau.) are the new spatial and time coordinates, respectively. By convention, the coordinates xi and tau appearing in the Region statement are uniformly spaced. In multidimensional problems, the space they define is logically rectangular. Next, we specify the equations defining the coordinate transformation to the "old" coordinates S and t, the coordinates in which the PDE's are defined. In general these equations are nonlinear relations, so that the uniform sampling of the logical space defined by the new coordinates is mapped into a non-uniform sampling of the "physical space" defined by the old coordinates. In multi-dimensional problems the logically rectangular regions of the new coordinates can be mapped into physical regions which may be of fairly arbitrary shape. Similarly, if the old time coordinate is given as a non-linear function of the new time coordinate, we can compress or expand time. Our spatial coordinate transformations here are time independent, i.e. stationary, but we can also define moving, or time-dependent coordinate systems. The System automatically transforms all the equations, interior, boundary, and initial, to the new coordinates. For example the Black-Scholes equation becomes: ##EQU12## where J(.xi.)=.differential.S(.xi./.differential..xi.) is the Jacobian. The transformation of this one dimensional equation is a fairly simple matter. In higher dimensionality however, such manipulations can become exceedingly complex, e.g. transforming the Navier-Stokes equations of fluid dynamics to three dimensional, time-dependent generalized coordinates. The System however, handles this readily. By default, equations are discretized in conservative form (the chained derivatives in the above equation are not expanded.) This is the preference of most experts, and is also appropriate when the coordinate transformation is not known analytically as it is here, but only in tabular form, as passed from a numerical grid generator, for example. In one space dimension, effective grid spacing in the original coordinate is proportional to J(.xi.) The transformation used here is ##EQU13## J(.xi.) and hence the grid spacing, is minimum near .xi..sub.0 3.4.2 Discrete Barriers To model a discretely sampled barrier, the System does not simply place the region boundary at the barrier and enforce the boundary condition of option value equal to 0 or to discounted rebate, if applicable. Because the boundary condition is enforced at every time step, this represents a continuously sampled barrier (even though time itself is discrete in our finite difference model). Instead, the System places the region boundary well beyond the barrier, S_{max}> Bin this case, let the option value "diffuse" over the barrier between sampling dates, and then enforce the knockout condition only at the discrete sampling dates. In the previous example, the time evolution template was not specified and hence, defaulted to a simple time loop with equal fixed time steps. The DiscreteEvents statement specifies a time evolution template in which the solution is integrated in time from the initial condition (option expiry) up to the first discrete event time. A discrete event, e.g., sampling of a barrier or path variable, or paying a discrete dividend, is executed by making the appropriate interpolation or assignment of option value and any hedging parameters. Then, integration proceeds to the next discrete event time. Each particular event (actually, series of events) has an associated array of event times. The generated code keeps a marker designating the next time .DELTA.t in each series and adjusts t, if necessary, to land exactly on the time of the next pending event, or if the last event has been executed, on the final time (the present in backward time.) There are currently four discrete event descriptors in addition to Barrier. They are Dividend, Asian, Lookback, and Path. Asian and Lookback are special cases of the more general Path. Dividend pays a discrete dividend. All four operate by making an interpolation in a spot price or path dimension based on path information accumulated at discrete intervals. Any number of descriptors may be specified in the DiscreteEvents command. In the Barrier descriptor, one gives the direction in which the barrier acts, here spot price S, a conditional representing the barrier crossing, and a list of functions or conditions to be applied when the barrier is crossed. Other arguments include files where the number of samples and the sampling times may be found. Optional arguments such as nsampl==nbar allow renaming a variable for the number of samples. The System then refers to this name in other parts of the specification file, and avoid name collisions. The direction specifications in Barrier seem superfluous in this single asset example, but they are needed in multi-asset problems where there may be different discrete barrier sampling times for each of several assets. 3.4.3 Numerical Results
TABLE 3.4.3
Comparison of finite difference and Monte Carlo
results for the value of an at-the-money European
call with a discretely sampled up-and-out barrier.
nbar is the number of discrete barrier samplings,
equally spaced. Other parameters are: TMax = 1.,
iMax = 5000, K = 100., nMax = 40000,
r = 0.05, sigma = 0.25, SMax = 200., B = 150.
DISCRETELY SAMPLED EUROPEAN KNOCK-OUT CALL
Monte Carlo Monte Carlo finite
Nbar 10,000 paths 100,000 paths difference
4 7.772 .+-. 0.118 7.922 .+-. 0.038 7.956
16 7.154 .+-. 0.112 7.294 .+-. 0.036 7.353
64 6.773 .+-. 0.108 6.904 .+-. 0.035 6.935
256 6.542 .+-. 0.105 6.667 .+-. 0.034 6.694
Table 3.4.3 compares results of the finite difference code with Monte Carlo results for the value of an at-the-money European call with a discrete knockout barrier as a function of the number of barrier samples. The finite difference grid is uniform and has 5000 grid points, reducing truncation error to insignificant levels. The spot and strike prices are S=K=100, the barrier is B=150. Other parameters are listed in the table caption. The finite difference and Monte Carlo results differ by about the Monte Carlo variance for all sampling intervals. However, the finite difference values are consistently higher. This is reasonable given the significant increase in the Monte Carlo values as the number of paths is increased. A further order of magnitude increase in paths would be expected to raise the Monte Carlo values, bringing them into closer agreement with the finite difference values. It is interesting that there is a large difference between even daily sampling, nbar=256 and the continuously sampled value of 6.442. FIGS. 7(a) and (b) show Results for a discretely sampled Europan knockout call, with: (a) present option value with a 200 point nonuniform grid, (b) effective grid spacing .DELTA. of the nonuniform grid, (c) present option value with nmax=1000 and nmax=10000, (d) error in the hedge parameter .DELTA. for 200 point uniform grid(solid curve) and 200 point nonuniform grid (dotted curve). Other parameters are: are Tmax=1., K=100, I=0.05, sigma=0.25, Smax=200, B=150, nbar=256. FIGS. 7(a) and (b) display finite difference results for the knockout call with nbar=256, i.e., approximate daily sampling. Just 200 grid points in spot price S have been used, but the coordinate transformation concentrates grid points near the barrier B. The present option value is displayed in FIGS. 7(a)(b) as a function of spot price. The value is small but finite for S.gtoreq.B, reflecting the possibility that the price will move below the barrier before the first sampling. The effective grid spacing, (2S/.differential..zeta.).DELTA..zeta. is displayed in FIGS. 7(a)) and (b). Using alpha=5, results in a ratio of about thirty for the largest to smallest spacing on the grid. The grid spacing near the barrier is about one seventh that of a uniform grid of the same dimension. While the Crank-Nicholson time differencing used in the examples is unconditionally stable, truncation error can still result in undesirable behavior if the time step is too large with respect to the characteristic grid diffusion time. When a coordinate transformation is used to reduce the grid spacing, the time step must also be reduced to avoid the problem. This is illustrated in FIGS. 7(a) and (b), which shows the present option value near the barrier for different numbers of time steps, nMax==1000 and nMax==10000. The larger value of nMax yields a smooth solution but with the smaller one, a spike appears near the barrier. This spike is not a symptom of numerical instability, but rather, that with the larger time step, the very sharp gradients that are periodically produced by setting the option value to zero for S.gtoreq.B on barrier sample dates, do not decay fast enough. Indeed, with Crank-Nicholson differencing, a Nyquist mode on a finite difference grid has an amplification factor of minus one as .DELTA.t.fwdarw..infin. and decays not at all. Although the spike remains localized near the barrier, it reduces diffusion through the barrier and results in an option value which is everywhere too large. The resulting error in value is much less than that for a uniform grid, but several times that for the properly time-resolved nonuniform grid calculation. The errors in the hedge parameters .DELTA. and .GAMMA. very near the barrier will obviously be worse than those in option value when the time step is too large. FIGS. 7(a) and (b) compare the error in the hedge parameter delta for uniform (solid curve) and nonuniform grids (dotted curve) of 200 points, properly resolved in time. (Here, error is defined relative to the results of 5000 point uniform grid results.) The effectiveness of the nonuniform grid is evident. The following Table 3.4.3 displays RMS errors and relative CPU times for 200 and 2000 point uniform grids and a 200 point nonuniform grid. The nonuniform grid reduces RMS errors by factors of between fourteen and twenty-five, relative to the 200 point uniform grid. The cost of the huge reduction in error is an increase in CPU time due to the larger number of time steps. However, this is still a factor of twenty less than the CPU time required to achieve similar accuracies with a uniform grid. Comparing the errors for the uniform grid results, it is evident that near the discretely sampled barrier, the nominal second-order accuracy of the default finite difference algorithm does not hold. 3.5 Stochastic Interest Rates Many interest rate derivative pricing problems can be posed as PDES, and are therefore suitable for synthesis. This example prices a cancelable swap under the Black-Karwinski (B-K) (see, F. Black and P. Karasinski, Bond and Option Prices When Short Rates Are Logonormal, Financial Analysts Journal, July-August, 1991, 52-59.) interest rate model, in which the interest process is given by dlnr=[v-lnr]dt+.sigma.dXx (Eqn. 30) where X is a Wiener process. P. Witmott, J. Dewynne, and S. Howison, Opti | ||||||