A Simulation of the Heston Model with Stochastic Volatility Using the Finite Difference Method

In this study, we investigated one of the most popular stochastic volatility pricing models, the Heston model, for European options. This paper deals with the implementation of a finite difference scheme to solve a two-dimensional partial differential equation form of the Heston model. We explain in detail the explicit scheme for the Heston model, especially on the boundaries. Some simple ideas to modify the treatment on the boundaries, which leads to a lower computational cost, are also stated. The paper also covers comparisons between the explicit solution and the semi-analytical solution.


Introduction
In recent years, financial markets containing derivatives have become more and more popular throughout the world. Derivatives were introduced in the Vietnamese equity market in 2017 and attracted a lot of attention. Options are derivative products that give the owner the right, but not the obligation, to buy or sell (depending on the type of option) an underlying asset at a stated price within a specific timeframe. One of the most concerning problems for traders is finding the fairest price for an option to incorporate into their strategies to maximize profits.
In the early 1970s, Fisher Black and Myron Scholes (Black-Scholes, 1973) derived an option pricing model that played an essential role in the development of modern derivatives finance. They used a partial differential equation to obtain values for European calls and put options on the stock. The model is computationally simple and quickly became one of the most popular models used to determine the value of an option.

Vietnam Journal of Agricultural Sciences
However, the assumptions made in the Black-Scholes model are rather restrictive. In particular, the Black-Scholes model assumes that the underlying asset price follows a geometric Brownian motion with a fixed volatility. The 1970s and 1980s observed many significant extensions of the model (Hull & White, 1987;Stein & Stein, 1991;Heston, 1993;Bates, 1996;Heston, 1997). Among the stochastic volatility extensions of the well-known Black-Scholes model, the Heston model, developed by Steven Heston in 1993(Heston, 1993 is the most widely used.
The Heston model assumes that the volatility   vt follows a stochastic process, which is a crucial improvement in comparison with the Black-Scholes model. According to Heston (1993), the price of an option ( , , ) U S v t must satisfy a two-dimensional convection-diffusion partial differential equation ( In the model above, there are several parameters that we need to determine and they will be introduced in more detail in the next section. Similar to the Black-Scholes case, these parameters can be calibrated using market data. This step involves big data collection and estimation methods. We also need to introduce the initial condition and boundary conditions for the Heston equation, which will be different depending on the options type. These conditions are also discussed below. In our work, we focus on European options. As in Heston (1993), the closed-form solution exists in an integral form, which is reviewed after the introduction of the model. For the other types of options, the analytical solutions are not expected to have explicit formulas.
Even though the closed-form solution exists, it appears in a proper integral containing complex functions. It leads to difficulties in computing the Heston solution. Therefore, a numerical method is usually used. Before reviewing the literature about PDE-based numerical methods, we want to mention that despite its slow process, the Monte Carlo method is also one of the most widely used numerical techniques for option pricing problems in general and for the Heston model in particular (Kahl & Jäckel, 2006;Lord et al., 2010). An advantage of the PDE representation is that this type of partial differential equation is well studied in the literature and many numerical methods have been developed to solve them. Since the Heston equation is a parabolic PDE, the standard finite difference schemes are quick approaches (Thomas, 1995;Hull, 2012;Rouah, 2013). Here, the explicit method is a popular approach due to its simplicity. However, the method requires a large number of time steps. Many authors instead use the implicit method or alternative implicit method (Douglas & Rachford, 1956;Peaceman & Rachford, 1955;Craig & Sneyd, 1988;Hundsdorfer, 2002;Hout, 2007;During & Miles, 2017) which do not limit the number of time steps but require at each time step the solution of large sets of equations. Other authors follow the ideas of the splitting method as a combination of operator splitting and iterative methods, and transform a 2dimenstional problem into quasi 1-dimensional ones (Safaei et al., 2018;Li & Huang, 2020). Alternatively, one can use the finite elements method (FEM) for the valuation of European options as in the papers of Winkler et al. (2001) or Chen et al. (2014). The FEM method is more complicated to implement but it has an advantage that it works for exotic options as well.
In our work, we expect to apply an effective method to determine the price of options in the developing Vietnamese derivative market. Since the Vietnamese derivative market is in its first steps of development, we chose the Heston model with a simple European option type to simulate. We first reviewed the numerical methods, which mostly depend on a Fourier transform, to deal with the implementation of the Heston solution (Rouah, 2013). In addition, we also followed one of the effective numerical methods to solve the Heston problem. Since the explicit finite difference method has an advantage of simplicity and it works well for parabolic PDEs if the stability condition is satisfied, we chose the explicit method and clearly showed how it works for the Heston model. Even though the approach is classical, we showed the schemes in every detail for this specific problem. A disadvantage of the explicit method is the requirement for the number of time steps. One of the main points in our work is the presentation of some ideas to modify the treatments on the boundaries specified in the remark of explicit finite difference of the Heston model, which leads to a more effective result in practice in the sense that the number of time steps is reduced, the size of the consideration domain is smaller, and at the same time, lower relative errors are obtained as shown in Table 6. The results are implemented in the Python programming language and we tested the accuracy of the method by comparisons with the Heston solution.
The structure of this paper is as follows. Firstly, we will review the Heston model and its closed-form solution given by Heston (1993). Secondly, we explain in detail the explicit scheme. Next, the algorithms to compute the Heston solution and explicit solution are discussed. We also include remarks about some ideas of the treatments on the boundaries. We show the comparisons of the explicit solution, explicit solution with the new treatments on the boundaries, and the Heston solution in the next section. We then provide further discussions about the stability and the rate of convergence. Finally, the last section covers the conclusion of all our work.

IoT-based hardware
As in the Black-Scholes model, the behavior of options on assets is assumed to follow the stochastic differential equation Sdt v t Sdz t   where S is the price of the underlying asset at time t ,  is the risk-free interest rate, () vt is the variance at time t , and 1 () zt is a geometric Brownian motion. In the Heston model, the volatility satisfies where 2 () zt is another Brownian motion,  is the reversion rate,  is the reversion level, and  is the volatility of the variance process. The correlation between the two Brownian motions is  , i.e, Following the method that works for the Black-Scholes model, with the help of Ito's lemma, one can demonstrate that the price of an option ( , , ) U S v t must satisfy the partial differential equation (Heston, 1993;Hull, 2012) (2.4) Here the term ( , , ) S v t  represents the price of volatility risk and is independent of the particular asset. In Heston (1993), the author also proposed that a European call option with strike price  and maturing at time T satisfies equation (2.4) with the boundary conditions and initial condition (2.5) ( , , ) max(0; ), (0, , ) 0, ( , , ) 1, Equations (2.4) and (2.5) are the governing equations for the Heston model. It is an extension of the Black-Scholes model in the sense that the behavior for the volatility process is assumed to be stochastic instead of being a constant. Note that the coefficients of the model need to satisfy the so-called Feller condition, 2 2 / ( ) 1    .
Throughout this work, we will consider the problems (2.4) and (2.5) to find the price of a European call option. Without a loss of generality, we can assume that 0.

 
A simulation of the Heston model with stochastic volatility using the finite difference method

544
Vietnam Journal of Agricultural Sciences

The closed-form solution of the Heston model
By analogy with the Black-Scholes model, in Heston (1993), the author derived a closedform solution to the problems (2.4) and (2.5) as follows. The solution has the form (2.6) () 12 ( , , ) , , Since the integrand in (2.7) decays rapidly, it is known that the integral is convergent and numerical computation of (2.6) is available. In our work, we will take the solution of the form (2.6) as the reference solution to show the accuracy of the numerical scheme.

A Finite Difference Scheme for the Heston Model
In this section, we will discuss finite difference schemes to solve (2.4) and (2.5) The problem is to find the approximations Now, we will explain in detail the treatments on the boundaries. The conditions on 0 S  and v  as (0, , ) 0, ( , , ) ,  At this step, the full problem of the Heston model for the European call option has been solved numerically by a finite difference explicit scheme. Our discussion above already covers the classical approach for this specific problem, moreover, we also clearly show the treatments on the boundaries.

Consistency
With the help of the Taylor series (Thomas, 1995;Hull, 2012), one can conclude that the explicit scheme is Thus, the scheme is consistent.

Stability
Following the standard Von Neumann analysis using Fourier techniques, the explicit scheme is conditionally stable. As in Sensen (2008), we take The condition for stability is very important to make sure that the scheme works. This is also one disadvantage among the many advantages of the explicit scheme.

Convergence
The convergence of this scheme is guaranteed by the Lax equivalence theorem. The basic idea is that, if a scheme is consistent, it is convergent when it is stable (Thomas, 1995). Thus, when the stability condition is satisfied, the explicit scheme is convergent.

Numerical Algorithms and Python Implementation
In our work, we will use the Heston solution of the form (2.6) as the reference solution. Since the closed-form solution (2.6) appears in an integral form, it also needs a numerical computation to find the analytical solution. Then, we will implement the explicit scheme for the Heston model.

Semi-analytical solution of the Heston model
The Fourier inversion transformation (2.7) is the main point in the numerical implementation of the Heston solution provided that the characteristic function is known. To fix the first problem, we have to truncate the semi-infinite integral to a finite one. In this algorithm, we chose a maximum level integrating the characteristic function j f , which is typically oscillating and causes the instability of the numerical scheme. Following a simple approach in Lord and Kahl (2007), we can replace j g in (2.7) with (4.1) namely, we just switch the minus and the plus signs in front of the j d . Lord and Kahl (2010) proved that with this modification, the numerical scheme to compute the options price is stable.
We will implement the Heston solution numerically as follows.
(1) Define a mesh on the space domain (2) Define a maximum level max U to compute the integral numerically.
(3) At each grid point ( , ) ij Sv of the mesh, compute The real coefficients ,, a b u as in (2.10), The complex coefficients , dg as (2.10) and (4.1), The coefficients , CD of the characteristic function as in (2.9), The integrand of (2.7), The approximation of 12 , PP using (2.7) and numerical integration with  varying from 0 to max U , and The approximation of We use the following parameters, as chosen in Sensen (2008), to implement the semianalytical solution of the Heston model. Figure 1 shows the result of the Heston solution with the parameters as in Table 1.
The implementation was produced in Python programming language. We have made use of the Numpy and Scipy libraries together with the cmath module which is very convenient when dealing with complex numbers as well as integration.

Explicit finite difference solution of the Heston model
The algorithm of the explicit scheme to find the approximation of the European call option via the Heston model is straight forward from the explicit scheme in the previous section . We also use the parameters as in Table 1.
Consider a discrete domain as shown in the previous section. Note that we need to choose a suitable time step t  such that the scheme is stable as in (3.13). In our work, we solve the problem in the domain We also take 5, 0.05 Sv with varying time step t  . Starting at time tT  , using the initial condition, the values at grid points are defined by Figure 2 shows the initial surface. Then, we compute the coefficients in (3.3), (3.7), and (3.11). At each time step, we will do as follows.
(1) For the points in the interior domain (without the boundary), the value   1 ,1 ,, The results of the explicit method are given in Figure 3 below.
Moreover, we also implement the explicit solution of the Heston model with the parameters as presented in Table 2.
The parameters in a real model can be determined using market data. This work requires a big data collection and an efficient Vietnam Journal of Agricultural Sciences   Table 1 Figure 2. The initial surface of the Heston model with parameters as in Table 1 estimation method. The problem of finding the parameters stands in other interesting fields and is not covered in our work. However, we expect to deal with this problem in the future as we apply the results to the Vietnamese derivative market. Here, the parameters in Table 2 are the ones estimated in Yang (2013) for the Heston model with the price the Google Inc. company recorded on April 6, 2013. We tried to solve the problem with a several values for expire time T and strike prices K . Figure 4 shows the results for 0.112328767 T  and 510 K  .

Remarks
(1) It turns out that the choice of t  is important to make the scheme stable. With the suggestion (3.13), we choose / 5000 tT    Table 1 Figure 4. The explicit result of the Heston model with the parameters as in Table 2 (2) Since the boundary conditions when v  and S  stated for large values of v and S , we need to take a bigger domain than the domain of consideration. In our work, we choose the space domain (3) A large number of grid points in the space domain enlarges the number of time steps, but it is expected that the domain of consideration should be as small as possible. One of the simplest ideas is that we can modify step (4)

 
, which says that the change of the price in the interior domain is "a bit" smaller than the change on the boundary.
All the numerical solutions are implemented in Python language with the help of the Numpy library, which provides a very powerful environment to work with arrays and vectors. Note that since the coefficients (3.3), (3.7), and (3.11) are independent of the time t , we compute them outside the time loop and employ a lower computation time.

Comparisons
We will take the Heston solution as the reference solution to compare. Let us take the first case to look at the errors in more detail. Figure 5 below shows the price surfaces given by both the Heston formula and the explicit scheme in the same coordinate system. The colored surface represents the explicit result, while the green mesh represents the Heston solution.
We also cut through the surfaces by the plane 1 v  to see the differences between two graphs as shown in Figure 6, where the higher green graph represents the explicit result and the lower purple graph represents the Heston solution. Figure 6 shows a small difference between the Heston solution and the explicit solution.
It is easy to see that the time step size t  strongly impacts the result. When t  is big, the scheme is unstable, and the error blows up as shown in Table 3. In contrast, when t  is small enough, the scheme seems stable with the error as illustrated in Table 4.
However, even though we expected that the error tends to 0 as 0 t  , it seems that the convergence is very slow and the error has small differences when t  varies from 1 1500 to 1 20000 . The relative error is given by Error As mentioned in the previous section, we need to use a bigger domain than the original one to cover the boundary condition when S  and v . To obtain a better error, we extend the domain of consideration to     ( , ) 0,300 0,1.2 Sv . The results are given in Table 5 below.
From the tables, we can see that the scheme has a smaller error when the size of the domain increases. We also can see that the error is getting smaller as time step 0 t  but the process is slow (first-order in time). On the other hand, the convergence is much faster when the space steps tend to 0 (second-order in space). However, if we take the space steps smaller, we need to take the time step extremely small due to the stability condition. It is a main problem of the explicit method in spite of its simplicity. As seen in . As can be seen in Table 6, we obtain a smaller error in comparison with the ones in Table 5. Moreover, the time steps t  are much smaller, starting from 1 2000 t  . Indeed, we can use the time steps as in Table 4. Therefore, our modification helps to reduce the size of the S  domain as well as the number of time steps, which gives a lower computational cost. The results are also better in the sense that the errors are reduced.
We also tried the Heston model with different parameters. Following Yang (2013), with the data of Google Inc., the parameters for the Heston model are chosen as in Table 2.   Table 1. Figure 6. The cut at 1 v  of the surfaces in Figure 5.

Conclusions
In our work, the explicit finite difference scheme with some simple modifications on the boundaries have been applied to the Heston pricing model with European options. Our results show that the explicit scheme with the modifications is simple, easy to apply, and gives a better approximation to the Heston solution in comparison with the classical method. However, the scheme still has the disadvantage of stability conditions. To overcome this disadvantage, some authors developed the alternative direction implicit (ADI) scheme, which is a kind of mixed explicit-implicit scheme. The ADI scheme is unconditionally stable. Nevertheless, it is more complicated in comparison with the explicit method. Furthermore, the ADI scheme must be used with care since a high dimensional matrix appears in a linear system of equations and an efficient solver is needed. In future work, we expect to apply this method in an efficient way to the Heston model with more complicated options types such as the American option and Barrier option. Besides, the parameter estimation and the application of this scheme in the Vietnamese derivative market is also a problem of interest.