A Lie Algebraic and Numerical Investigation of the Black-Scholes Equation with Heston Volatility Model

This work deals with an extension of the Black-Scholes model for rating options with the Heston volatility model. A Lie-algebraic analysis of this equation is applied to reduce its order and compute some of its solutions. As a result of this method, a five-parameter family of solutions is obtained. Though, these solutions do not match the terminal and boundary conditions, they can be used for the validation of numerical schemes.


Introduction
Black and Scholes [1] assumed a financial market, where a risk free bond with constant interest rate r, an asset with price S that is modelled by a geometric Brownian motion, and call and put options related to this asset can be traded. With the assumption of an arbitrage free market and in the framework of It o 's stochastic differential equations, it is possible to derive the well-known Black-Scholes partial differential equation for the fair price of an option V. As the assumptions of this first modelling attempt are in practice too restrictive, several extensions of this model were proposed. One possible direction is to discard the assumption of a constant volatility for the geometric Brownian motion of the asset price and to assume that it is itself a random variable governed by a stochastic differential equation [2]. The resulting stochastic differential problem is given by (1) = , where we suppose that (1) (2) , t t W W are two stochastically independent Wiener processes. The model constants α, m and L are supposed to be positive. The drift term of (1b) is built in a way such that the average of v t tends to have approximately the value m. In particular, if L is zero then v t is deterministic and converges exponentially to m as t tends to infinity. In this case, the option price behaves according to the solution of the usual Black-Scholes equation with a constant volatility = m σ . Based on (1), it is possible to derive a partial differential equation for the price of an option V [2] as it has been done for the model with constant volatility in studies of Gunther and Jungel [3]. In contrast to the standard Black-Scholes equation the PDE that arises with Heston's volatility model involves one more argument representing the current volatility of the market. The resulting two-dimensional Black-Scholes equation is as follows.
In application, this equation is augmented with the following terminal and boundary conditions In this setting, a denotes whether a call (a = 1) or a put (a = −1) option is considered. T represents the time when one is allowed to buy or sell a share of an asset for the prescribed price K, whereas λ is the parameter that models the price of volatility risk [2]. In (2), x represents the asset price S and y denotes the current volatility v. As the asset price and the current volatility are always positive, we are searching for a solution of the above differential problem (2)-(3) in the domain In the following, analytical and numerical solutions of (2) together with the boundary and terminal conditions (3) are sought. In particular, a quick review of Lie symmetries of partial differential equations is given in Section 2. In Section 3, this method is applied to the 2-dimensional Black-Scholes equation (2) and we derive a five-parameter family of analytical solutions. In Section 4, convergence properties of the Chang-Cooper discretization are tested with the given analytical solutions. A section of conclusion completes the exposition of our work.

Lie Theoretical Analysis of Differential Equations
In this section, we illustrate how Lie symmetries can be used to determine analytical solutions of partial differential equations. Applications of this method can be found in literature of Bordag [4] and Naicker V, Andriopoulos K, Leach [5]. Our review is based on the book of Stephani [6].
Many partial differential equations for a function u that is dependent on n variables x i ( = 1, , i n  ) can be written as follows with an analytic function H, where y k denotes subsequently the independent variables x i , the dependent variable u, and its derivatives where the coefficients 1 i i s η  are given by Here, the total differentiation operator A Lie symmetry of a PDE is defined as a group of transformations of the independent and dependent variables such that set of solutions is invariant under these transformations. From the fact that the image of a solution satisfies the PDE, i.e. ( ( )) = 0 k H f y ε for all ε, it can be shown that holds everywhere on the solution manifold H(y k ) = 0, where X is the prolonged symbol of the transformation group. For a given Lie symmetry, we seek to find its canonical variables w k as the corresponding symbol is  and (9) then reads as . Hence, the resulting PDE written in the new variable w k is independent of w 1 and therefore involves one independent parameter less. Computing solutions of the transformed PDE, which should be easier as less independent variables are involved, and reversing the transformation, provides solutions of the original PDE.

Lie Analysis of the Two-dimensional Black-Scholes Equation
In this section, the Lie method is used to find solutions to the 2-dimensional Black-Scholes equation (2) written as H = 0, where H is defined as follows .
Then, we first apply the prolonged symbol to the function H and then evaluate the resulting function X(H) on the solution manifold H = 0. The resulting expression yields zero, whenever X is the generator of a Lie symmetry. The exact expressions for the prolonged coefficients ϕ 1 , ϕ 2 , ϕ 3 , ϕ 11 and ϕ 22 according to (8) are given by where sub-indices x, y, and V of ξ, γ, τ and ϕ denote partial derivatives with respect to the given variables. The equation H = 0 is solved for V t and inserted into X(H) = 0. Afterwards this single equation splits up into the determining equations, since the derivative variables (V x , V y , V xx , …) are linearly independent. Among the resulting equations, there are simple ones as . Hence, we solve the following remaining system of partial differential equations Notice that equations (IV,V and VI) are similar to the original PDE we are trying to solve. Inserting the special form of ϕ into (VI), it splits up into the following two equations due to the fact that φ, β, and τ t are independent of V. The function β is independent of the other functions and equations. Furthermore, it must be a solution of the PDE (VI'') whose Lie symmetries we are looking for. So the transformation V V εβ is a Lie symmetry, mapping solutions onto solution. As the Lie symmetry corresponding to the coefficient function β is as difficult to find as solving the PDE directly, it is not significant for our purpose and we do not take (VI'') into account any more.
Instead, let us focus on (II). Differentiating it two times with respect to y yields 3 , whose solution is given by . Inserting this expression in (II) gives γ = τ t y. Hence, γ is independent of x and since (I) holds, the function ξ is independent of y, i.e. = = 0 The same idea works with (III), whose second derivative with respect to x is The general solution to this ordinary differential equation With this knowledge, Equations (III) and (IV) simplify to Notice that the functions τ and D are independent of y.
Using previous results in (VI'), we obtain As the coefficient of y must equal zero, τ t = 0 holds and consequently B t is equal to zero. Hence, the most general solution of the determining equations is where c 1 , c 2 and c 3 are real constants and β is a solution of the Black-Scholes equation.
Hence, the only Lie symmetries that (2) admits have the following infinitesimal generators Apart from the last symmetry this is a three dimensional, solvable Lie algebra, i.e. the commutator of two arbitrary symmetries X and Y equals zero.
Next, we determine the canonical variables for the Lie symmetry with fixed constants c 1 , c 2 and c 3 . Therefore, we search for three functionally independent invariants ( , , , ) v x y t V , u x y t V that satisfy the following equation The evaluation of (7) with the choice of the following new variables , and log( ) shows that the Lie symmetry with respect to the new variables has the desired normal symbol ∂ ∂ and that v, w, u, and s are the canonical variables.
In order to rewrite (2) in the new variables, we differentiate with respect to x, y and t and obtain x v xx vv Hence, the 2-dimensional Black-Scholes equation in the new variables is given by and by setting 1 3 2 = = c c rc we cancel out terms with u and u v . Therefore, the reduced Black-Scholes equation is as follows In order to find solutions to (10), we assume ( , ) = ( ) ( ) Since the left-hand side of the equation depends only on w and the right hand side only on v, both sides must be equal to a constant C. Hence, we obtain two decoupled ordinary differential equations The general solution of equation (11) is given by Regarding the second ordinary differential equation (12), we transform it into Kummer's equation To summarize, we obtain the following five parameter family of solutions to the two-dimensional Black-Scholes equation   So if V would be a solution to the two-dimensional Black-Scholes equation subject to the given boundary conditions and if it corresponds to a solution u of the reduced equation, then u must satisfy the boundary conditions above. Note that, while (i) imply (ii) and (iii), (i) is not consistent with (iv). As (i) determines u and therefore u is independent of w, (iv) cannot be satisfied as

Numerical Solution of the Two-dimensional Black-Scholes Equation
This section deals with a numerical scheme to calculate an approximation to the solution of the Black-Scholes equation (2). We work with the proposal of Chang-Cooper scheme [9] and analyzed in studies of Mohammadi and Borz [10]. This disretization scheme is often used for Fokker-Planck equations, as its solutions are probability density functions and therefore are non-negative and their integral over its domain equals 1. These two properties are preserved by the Chang-Cooper (CC) difference scheme. In the case of the Black-Scholes equation the solution is also non-negative, as it models the price of an option, which must be non-negative. Hence, the choice of the CC scheme guarantees that the numerical solution will be non-negative. In order to apply the Chang-Cooper discretization scheme the two dimensional Black-Scholes equation (2) we obtain the following PDE We can write (15)   domain. This is the case when the condition is satisfied.
The transformed Black-Scholes equation (16) must be solved subject to the following transformed initial and boundary conditions:  with its coefficients are added to the right hand side of the linear system of equations, as they are known.
• The boundary condition for x → ∞ is given only in terms of the derivative of V with respect to x. Therefore the value • The derivative term is known and can be put to the right-hand side of the equation.
The values of the function on the boundary y = 0 are not given. Fortunately, as y goes to zero, B y goes to In the following, the numerical scheme is applied to the function type C = 0 in (14). After the variable transformation, that is used to write (2) in flux form, the test function becomes  The plot data is shown in Table 1 where N, M and Q is the number of grid points in the x-, y-and t-dimension, respectively. A small time-step size is used in order to have a small error for the time discretization and to investigate the dependence of the error on the spatial-grid size.
The numerical experiments show that the discretization that is used provides only first-order convergence, i.e., doubling the grid point number in each spatial dimension and therefore halving the grid size h results in an error that is half as big as before. Notice that secondorder convergence is proven in literature of Mohammadi and Borz [10] with the assumption of zero boundary conditions. In order to validate this theoretical result the same procedure is done with a Gaussian bell function that is almost zero on the boundaries: In contrast to the test function f, this function Φ does not satisfy the partial differential equation (16 is added to the right hand side of the linear system of equations according to ( ) n  m  i  i  m  m  m  I  I  I  I  I  I Numerical experiments are performed with the grid sizes according to Table 2.  Figure 2 we can see that second-order convergence is obtained.
In the following, we use the Chang-Cooper numerical scheme to calculate a numerical solution. Here, this numerical solution in the case of a call option is compared to the solution of the Black-Scholes equation, where the volatility is assumed to be constant. It is given by where Φ is the cumulative distribution function of a normally distributed random variable with mean 0 and variance 1. This function is given by The model constant m represents the square of the average volatility and the stochastic process tends to this value. Hence, if one starts the process with the value y = 0.5, the stochastic process for the volatility t ν is likely to be almost constant to 0.5. As you can see in the lower left plot of Figure 3, the calculated price is nearly the same with both models. In contrast, regarding the case of a currently volatility lower than 0.5 the price of the option calculated with the extended Black-Scholes equation is higher than that of the model with constant volatility, because it takes into account that the volatility will rise. This In addition, the numerical solution satisfies the so-called Put-Call-Parity. The price of a call option C and the price of a put option P subject to the same asset with price x, that have the strike price K and the expiry date T in common, are related by the following formula [3] ( ) = .
r T t x P C Ke − − + − We compute also the price for the put option and observe the absolute deviation for the Put-Call-Parity formula that we average along the y-dimension. Figure 4 shows the result depending on x and t. Apart from small x values the error is in the range of the numerical error of the Chang-Cooper-Scheme. The drastic increase of the error for x → 0 is due to the fact, that the boundary condition for = To conclude, it is evident why there is such a great error for small x, and moreover it is not relevant as x gets never so small in applications.

Conclusion
The aim of this work was to solve the partial differential Black-Scholes equation with Heston volatility model. Therefore, an analytical technique due to Sophus Lie that can be use to reduce the number of independent variables of a partial differential equation was presented and applied to the Black-Scholes equation. A five-parameter family of solutions was found. These functions do not satisfy the boundary conditions of the option price problem and henceforth numerical schemes are necessary to obtain approximate solutions. In the last part of this work the Chang-Cooper discretization scheme was used to calculate the option price function numerically. Its convergence was tested with an exact solution of the PDE, which was found by the Lie theoretical analysis. Finally, the numerical scheme was applied to compute the price of an option and good result were obtained in accordance with economic reasoning.