This text is a piece on approximation of various stochastic models, with particular emphasis on schemes that avoid impossible solutions, for example negative volatilities.

** Comments and feedback welcome **

Some schemes we discuss are:

1 Euler-Maruyama method

2 Milstein methods (various forms)

3 Fitted Euler methods (Mirani/Duffy)

4 Predictor-Corrector methods

There are a number of popular finite difference schemes that are used the approximate the solution of SDEs. The most famous method is the Euler-Maruyama scheme and it has been applied in finance applications. It is applicable to a wide range of one-factor and n-factor models, including geometric, mean-reverting, square-root and stochastic volatility models. It is robust and easy to program but has a low order of convergence and the results that it produces are wrong in a number of cases. In particular, the methods can non-positive values and this is clearly impossible in financial applications because the quantities being modeling are typically stock prices, volatility or forward rates, for example. For this reason, a number of methods have been proposed that fix some of these inherent limitations of the Euler-Maruyama (EM) method. The most famous one is the Milstein method and it is in essence EM but with some added diffusion to compensate for the errors with the EM method. Again, the standard Milstein method can produce non-negative solutions and we have further noticed (from numerical experiments) that it is not more accurate than EM. This is well-known and a number of generalisations have been propsosed, for example the Balanced Implicit Method (BIM) and the Balanced Milstein Method (BMM). We now discuss a variation of the EM whose motivation is based on the so-called exponentially-fitted difference schemes that were developed by de Allen and Southwell in the 1950?s and elaborated by researchers such as Il?in, Emel?janov, Miller and many others (including some of my own work on convection-diffusion equations) thereafter. The idea is very simple and the method is based on fitting the drift and diffusion parameters for a general SDE with non-constant and even non-linear terms based on the fact that we know the closed-form solution for a one-factor SDE with constant coefficients. The method works quite well and in particular it gives excellent results for square-root processes. It is also possible to define a fitted Crank-Nicolson scheme and this is particularly good for mean-revertin*g models. We have not examined fitted methods for n-factor models. In general, BIM, BMM and fitted Euler schemes compensate for the fact that the Wiener increments can have both negative and positive values at each time level. They embody the structural properties of the SDE that they approximate. Finally, the Predictor-Corrector seems to be the most robust and applicable method for SDEs based on our numerical results and mathematical analysis to date. Actually the Predictor-Corrector method is explicit - we do an explicit Euler in the predictor stage of the scheme to get an intermediate value and then use this in a modified Trapezoidal rule. We have used it on a number of 1 and 2 factor SDEs and the results are very robust in general. This conclusion is based on numerical experimentation and it now a good moment to examine the scheme based on mathematical analysis. To this end, we wish to predict or have a formula that will allow us to determine the mesh-size (sometimes called h, k or deltaT) that ensures that the scheme has a unique, non-negative convergent solution. We apply the Banach fixed point theorem (BFPT): let (X, d) be a complete metric space and T: X -> X be a contraction operator. Then there is a solution of x = T(x) in X. For the current problem, we can write the Predictor-Corrector method in this for X =F (X) at each time level. The mapping F will be a contraction if we choose the mesh size k small enough and this is usually calculating the norm of the Jacobian matrix (for example, the spectral radius). I still have to some experiment to determine which norm is optimal in some sense.

The main features of the Predictor-Corrector method are:

. Functionality: it is applicable to a wide range of one-factor and n-factor linear and non-linear SDE. Good accuracy.

. Reliable: we have found it to be very robust in general. It is used in many areas of Numerical Analysis, for example PDE splitting methods

. Usable: it is relatively easy to program, especially if we have a generic framework; we only need to define the vector-valued functions for the drift and diffusion terms