We know that the solution for the propagator is given by
(4.57). Hence, the option price is given by
The error of the Monte Carlo method goes as as a consequence of the central limit theorem as long as there is no correlation between the configurations produced. (Though in general this condition is difficult to satisfy, we shall see later that we can easily satisfy it for this case). While this error may not look very impressive, it is often the best that can be managed for which have a large number of dimensions. For the present problem, we have a very large number of dimensions (in fact the exact problem has an infinite number of dimensions) and the Monte-Carlo method is the most practical one available for it.
We now have to produce configurations with the probability
distribution . While looks rather complicated, it has a
simple interpretation. It is the probability distribution for a
discretised random walk performed by . To see this, let us first
use Ito's lemma to find the process followed by . We find
The above result can also be derived by considering (4.6) as Brownian motion on a one-dimensional Riemmanian manifold with metric tensor .
Hence, in this simulation, we will use Euler's method to find the volatility paths since these paths are generated with the requisite probability distribution.
Hence, the algorithm to find a Monte Carlo estimate of the propagator is as follows
The propagator must finally be integrated over the final wave function to obtain the option price. The accuracy of this numerical quadrature depends on the spacing between successive values of . This means that we have to find the propagator for several values of to obtain reasonable accuracy which is computationally very expensive. We found it better to find the propagator using the above Monte Carlo method for only about 100 equally spaced values of over the range of the quadrature and using cubic splines to interpolate it at the other quadrature points. This produces excellent results (as we shall see later) as the propagator seems to be an extremely smooth function of .
Hence, the revised algorithm is of the form
For the special case of functions which can be expressed as piecewise polynomial functions, the integration can be carried out analytically and we can then just generate the paths and find the expectation over the generalised Black-Scholes prices. While this method has the advantages of being elegant, easy to program and completely eliminating the error due to the integration, it has the disadvantage that we have to carry out the integration first before we can write the program which will be specific to the given payoff. We have implemented this approach for call options (the code is given in appendix C) but found no significant improvement in run time (it runs about 1% faster).
The astute reader will note that there is one problem with the above algorithm in that it is unstable as since in this case . In the case of our simulations this is not a problem as our step size was always large enough to avoid this problem. If this is a problem in any simulation, this can be easily handled as the propagator can always be written as a functional of the paths as given in (4.64). We can then generate paths for or using (4.6), store the functionals (we will now need 5 terms (which is the same as above), , , and ). We can then proceed to find the propagator as above (of coure, now writing as a function of ) and integrate over the final payoff. For the special case of the European call option, one can make this method slightly more efficient by averaging over the generalised Black-Scholes prices (4.65). We did not use this method since it takes longer to compute the propagator using five terms and because the memory requirements are higher.
The normal way of doing Monte-Carlo simulations to find option prices
with stochastic volatility involves directly simulating the process
by discretising it using the Euler method to
where the standard normally distributed random variables and with correlation are generated by first generating two uncorrelated standard normally distributed random variables and and using and . The final values of are stored and the option price with strike price is estimated by considering . This is the algorithm used along with the control variate method in Johnson and Shanno (we henceforth call the above algorithm the standard algorithm).
Before we compare the efficiency of our algorithm with the standard one, we look at the possible sources of error and how the error due to each source scales with run time for both. The major source of error for both algorithms is the Monte Carlo error which goes as and which goes as the square root of the run time. Another common source of error in both the algorithms is the discretization of the process for or . The error in is then of the order of where is the time step used. However, the effect of this error on the option price is virtually impossible to estimate and we verify that the number of time steps is adequate empirically by comparing two simulations with different numbers of time steps. The standard algorithm has a further error due to the discretization of the process. The error is again of the order of . The effect of this error on the option price is of the same order as the payoff of the option is piecewise linear with respect to . Hence, this error also goes inversely as the square root of the run time. Our algorithm has other errors due to the finiteness of the integration over , the interpolation error and the quadrature error which goes as . The error due to the finiteness of the quadrature domain can be made completely negligible with minimal effort as the propagator goes as where . Hence, for large enough , this error goes as where is the run time. This is a truly negligible error. The interpolation error is difficult to estimate but we empirically show that it is very small for interpolation over 100 points as the propagator is very smooth. The quadrature error goes as (Simpson's rule). Since most of the computer time is spent on the Monte Carlo simulation, we can also make this error very small with only a small increase in run time. Hence, we see that one of the main advantages of our algorithm is that it has a significantly lower error for the same number of configurations (since it has no error due to the simulation since these degrees of freedom have been integrated out).
When generating a single set of option prices with the same error tolerance, our algorithm is about 30 times faster. However, our algorirhm has an important advantage in that and are independent of . Hence, when we generate sets of data with all parameters except fixed, we only need to calculate the new propagator using the terms and integrate over the final payoff. This effectively results in an increase in efficiency of a factor of six to seven when we calculate the prices for 10 different values of . Hence, when generating data with several values of , our algorithm is about two hundred times faster.
For 50,000 configurations and 128 time steps, the time taken by our program on a Silicon Graphics with a MIPS R4400 processor running at 250MHz is 96.8s and the standard error for an at the money 90 day option when , , , and is $0.006. For a similar simulation using exactly the same parameters and initial values, the standard algorithm takes 91.0s and has a standard error of $0.038 for the at the money option. Thus, for a similar run time, the standard error using our algorithm is about a factor of 6 smaller. Hence, to achieve a similar accuracy, the standard algorithm takes about 36 times longer. (The program involving the generalised Black-Scholes prices runs in 95s for the same values of the parameters.)
The control variate method employed in Johnson and Shanno reduces the variance in the standard algorithm by a factor of about 20 making it approximately that many times faster. In that case, we see that our algorithm is only slightly better when calculating one set of option prices, but is still significantly when simulating several sets for different .
One straightforward test that we applied was to check that our program gave the same results for zero correlation as the straightforward solution presented in chapter 5. Since the volatility paths must also be generated for the straightforward solution we used the same random number seed, time steps and number of configurations. This provides a rigorous test that the intermediate errors in our program (due to the finiteness of the integration domain and the numerical quadrature) are well under control. We present two of these tests, one with the same random number seed and one without in figure 5.1. We see that most of the error is due to the Monte Calo simulation rather than the integration over the final payoff. When we consider the fact that the Monte Carlo error is significantly higher for non-zero , this shows that the interpolation and quadrature error are very small.
For the non-zero correlation case, we have also compared our simulation results with those in Heston and Johnson and Shanno and have always found agreement between our results and theirs within two standard errors (the 95% confidence interval). We have also compared the results of our program with the results from the ``standard algorithm'' presented in the previous chapter. Hence, we conclude that our algorithm is essentially correct.