Merger J* and Borzi A
Research Assistant, Department of Mathematics, University of Würzburg, Chair for Scientific Computing Emil-Fischer-Straße 30, Room 02.013, 97074 Würzburg, Germany
Received date: December 19, 2015; Accepted date: January 25, 2016; Published date: January 27, 2016
Citation: Merger J, Borzi A (2016) A Lie Algebraic and Numerical Investigation of the Black-Scholes Equation with Heston Volatility Model. J Generalized Lie Theory Appl S2:006. doi: 10.4172/2469-9837.1000S2-006
Copyright: © 2016 Merger J, et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Visit for more related articles at Journal of Generalized Lie Theory and Applications
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.
Lie algebra; Black-Scholes equation; Differential equations; Lie symmetries; Diffeomorphisms
Black and Scholes  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 ’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 . The resulting stochastic differential problem is given by
where we suppose that 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 vt tends to have approximately the value m. In particular, if L is zero then vt 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
Based on (1), it is possible to derive a partial differential equation for the price of an option V  as it has been done for the model with constant volatility in studies of Gunther and Jungel . 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 . 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.
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  and Naicker V, Andriopoulos K, Leach . Our review is based on the book of Stephani .
Many partial differential equations for a function u that is dependent on n variables xi () can be written as follows
with an analytic function H, where yk denotes subsequently the independent variables xi, the dependent variable u, and its derivatives Equation (4) defines a manifold in some multi dimensional Euclidean space and its solutions are sub-manifold. Diffeomorphisms of Rn can be used to permute the set of solution and find solutions with special properties. Therefore, we observe that a one parameter group of diffeomorphisms can be completely determined by the first order differential operator
which is called infinitesimal generator or symbol. In fact, the group action on a point can be computed by solving the initial value problem (6)
Moreover, the coefficients of the symbol with respect to a change of variables y = T(x) with a transformation can be computed as follows
where X(Ti) denotes the application of the first order differential operator X on the function Ti. It can be shown that there always exists a set of canonical variables in which the symbol has the normal form
As the symbol X acts only on the independent variables, it is prolonged to act in a higher dimensional space including also the dependent variable and its derivatives up to the order of the partial differential equation. The prolonged infinitesimal generator is defined as follows
where the coefficients are given by
Here, the total differentiation operator is used.
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. for all ε, it can be shown that
X(H(yk)) = 0, (9)
holds everywhere on the solution manifold H(yk) = 0, where X is the prolonged symbol of the transformation group. For a given Lie symmetry, we seek to find its canonical variables wk as the corresponding symbol is and (9) then reads as Hence, the resulting PDE written in the new variable wk is independent of w1 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.
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
We assume that this equation admits a Lie symmetry with the infinitesimal generator
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 Vt and inserted into X(H) = 0. Afterwards this single equation splits up into the determining equations, since the derivative variables (Vx, Vy, Vxx, …) are linearly independent. Among the resulting equations, there are simple ones as Hence, we solve the following remaining system of partial differential equations
for the functions ,, and . 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 with its infinitesimal generator 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 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.
The same idea works with (III), whose second derivative with respect to x is The general solution to this ordinary differential equation (ODE) is Notice that (III) can only be satisfied by With this knowledge, Equations (III) and (IV) simplify to
From (IV') one can directly derive as φy is not dependent on y. Differentiating (III') with respect to x gives
Notice that the functions τ and D are independent of y. Hence
and hold, i.e. (III') becomes which means that andprevious results in (VI'), we obtain
As the coefficient of y must equal zero, τt = 0 holds and consequently Bt is equal to zero. Hence, the most general solution of the determining equations is
where c1, c2 and c3 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 c1, c2 and c3. Therefore, we search for three functionally independent invariants , , and that satisfy the following equation
The evaluation of (7) with the choice of the following new variables
, and 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
Hence, the 2-dimensional Black-Scholes equation in the new variables is given by
and by setting c1 = c2 = rc we cancel out terms with u and uv. Therefore, the reduced Black-Scholes equation is as follows
In order to find solutions to (10), we assume and obtain
which is equivalent to
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
Regarding the second ordinary differential equation (12), we transform it into Kummer’s equation
by defining where
and The general solution of Kummer’s equation is given by
where M and U are Kummer’s functions of the first and second kind, respectively. For further details, . In case of C = 0, (12) is given by
and we directly see that the first derivative of Ψ is a multiple of
and hence a general solution in the interval (0, ∞) is given by
Having found solutions u for (10), we obtain solutions V of (2) by applying the reverse variable transformation as follows
To summarize, we obtain the following five parameter family of solutions to the two-dimensional Black-Scholes equation
These functions do not satisfy the boundary conditions (3) given in Section 1. In order to check this we write down the boundary conditions in the new variables v, w, and u. They are given by
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 which equals neither v nor
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  and analyzed in studies of Mohammadi and Borz . 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) must be written in flux form. This is not possible, as the coefficient of V is −r and not However, introducing the following new variables ,, and and computing the derivatives of V with respect to the new variables as follows
we obtain the following PDE
We can write (15) in flux form as follows
At this point we would like to stress three important properties of the flux functions. To begin with, they are all independent of the time variable t. Hence, the left-hand side of the resulting linear system of equations is the same for each time iteration and the corresponding matrix must be computed only once. Moreover, both B x and B y are linear functions and therefore Lipschitz continuous with the constants
and Finally, these functions must be positive in our domain. This is the case when the condition
The transformed Black-Scholes equation (16) must be solved subject to the following transformed initial and boundary conditions:
There are several problems that arise during implementation:
• The domain of the problem (16) subject to (18) is and therefore unbounded in the space dimensions. Moreover the boundary conditions are given as a limit. For numerical purpose the domain was limited to and it was assumed that the function attain the limit values already at the finite boundaries.
• When corresponds to a point outside of the domain the values 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 is approximated by
• 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, By goes to and tends to zero. Assuming ,tends to zero, as goes to infinity. Hence, the coefficient of is zero and this function values need not to be known for the calculations .
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
Note that this function has a singularity in y = 0, if and only if .As this infinite value might arise problems while numerical calculations, the set of parameters is chosen such that holds. In particular, the test function was calculated in the domain with the following parameters
Unfortunately, it is not possible to chose a set of parameters such that the test function has no singularity at y = 0 and additionally there exists a domain where all flux functions are positive. That is because imply and therefore the necessary condition for (17) is not fulfilled. Consequently, there is no proof in this case, that the numerical solution is positive. Nevertheless, the convergence order can be observed. Figure 1 shows the difference of the numerical solution to the exact test function in terms of the norm
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.
|hx = hy||1/25||1/50||1/75||1/100|
Table 1: Numerical error for different grid sizes.
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  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). Hence the deviation to the PDE
is added to the right hand side of the linear system of equations according to
Numerical experiments are performed with the grid sizes according to Table 2.
|hx = hy||1/25||1/50||1/75||1/100|
Table 2: Numerical error for different grid sizes.
Φ was approximated in the domain with the parameters being
Here, the flux functions are positive, as With the logarithmic plot of the error in 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 spatial domain of discretization is After reversing the variable transformation the option price can be evaluated for and The following parameters are used
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 is likely to be almost constant to 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 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 effect can be seen in upper left and right plots. Finally, in the lower right plot the simpler model overestimates the price, when the initial volatility is higher than , due to the fact that it is likely to fall.
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 
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 for the Black-Scholes equation in flux form is applied at the finite value .This corresponds to as the transformation was Consequently, the numerical solution for the put option takes the value at whereas the correct value is as C tends to zero as x goes to zero and therefore the price for the put option is 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.
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.
This work was supported in part by project “Multi - ITN STRIKE - Novel Methods in Computational Finance” by EU Grant Agreement Nr. 304617 and by BMBF project 05M2013 ‘ROENOBIO: Robust energy optimization of fermentation processes for the production of biogas and wine.