Research Article  Open Access
Cubic Spline Method for a Generalized BlackScholes Equation
Abstract
We develop a numerical method based on cubic polynomial spline approximations to solve a a generalized BlackScholes equation. We apply the implicit Euler method for the time discretization and a cubic polynomial spline method for the spatial discretization. We show that the matrix associated with the discrete operator is an Mmatrix, which ensures that the scheme is maximumnorm stable. It is proved that the scheme is secondorder convergent with respect to the spatial variable. Numerical examples demonstrate the stability, convergence, and robustness of the scheme.
1. Introduction
An option is a tradable financial contract whose value depends on the value of an underlying asset. The buyer of the contract obtains the right, but not the obligation, to buy or to sell an asset at a specified price on or before a maturity date. A call option provides the right to buy the underlying asset for a certain price, whereas a put option confers the right to sell the underlying asset for a certain price. A European option can only be exercised at the maturity date, while an American option can be exercised at any time prior to its maturity date. Black and Scholes [1] showed that the value of a European option is governed by a secondorder parabolic differential equation with respect to the underlying asset price and time, which is known as the BlackScholes equation. The value of an American option is determined by a linear complementarity problem or as a free boundary value problem involving the BlackScholes differential operator [2]. It is often necessary to use numerical methods to solve these partial differential equations, as analytic solutions are not generally available.
There are several numerical methods for the valuation of European and American options. The first numerical method for option pricing was the lattice method proposed in Cox et al. [3] and was improved in Hull and White [4], which is equivalent to an explicit timestepping scheme. Since the BlackScholes equation with constant or spaceindependent parameters can be transformed into a diffusion equation, the finite difference methods applied to constantcoefficient heat equations have also been developed (see, e.g., Schwartz [5], Courtadon [6], Wilmott et al. [2], and Rogers and Talay [7]) for pricing options. Vázquez [8] presented an upwind numerical approach for the BlackScholes equation. Cen and Le [9–11] presented stable finite difference schemes on a piecewise uniform mesh for pricing European and American options. Wang [12] and Angermann and Wang [13] proposed a fitted finite volume method for the discretization of the BlackScholes equation. Rambeerich et al. [14] applied the exponential time integration scheme to price options. Other methods, such as meshless approach [15, 16], elementfree kpRitz method [17, 18], and elementfree Galerkin method [19], also can be used to solve the generalized BlackScholes equation.
The spline collocation methods are useful methods for solving partial differential equations. Spline solutions have their own advantages. For example, once the solution has been computed, the information between mesh points is available. Numerical methods based on spline collocation methods also have been used to solve option pricing problems. Christara et al. [20] proposed a quadratic spline collocation method to the American option pricing problems. Holtz and Kunoth [21] developed a Bsplinebased monotone multigrid method for the valuation of American options. Khabir and Patidar [22] applied a Bspline collocation method to solve the heat equation which is obtained from the BlackScholes equation by an Euler transformation. Kadalbajoo et al. [23, 24] used cubic Bspline collocation methods for the BlackScholes equation.
In this paper, we present a numerical method based on cubic polynomial spline to solve a generalized BlackScholes equation. We combine the implicit Euler method for discretizing the time variable with the cubic polynomial spline scheme for discretizing the spatial variable. The matrix associated with the discrete operator is an Mmatrix, which ensures that the scheme is maximumnorm stable. We will show that the scheme is secondorder convergent with respect to the spatial variable. Numerical results support the theoretical results.
The rest of the paper is organized as follows. In the next section, a generalized BlackScholes equation is introduced. The temporal semidiscretization is described in Section 3. The spatial discretization is constructed in Section 4. The fully discrete scheme is presented in Section 5. Finally, numerical experiments are provided to support these theoretical results in Section 6.
2. The Continuous Problem
We consider the following generalized BlackScholes equation: equipped with the terminal and boundary conditions Here, is the value of European call option at asset price and at time , is the exercise price, is the maturity, is the riskfree interest rate, and is the volatility function of underlying asset. When and are constant functions, it becomes the classical BlackScholes model. The existence and uniqueness of a classical solution of (1)(2) are well known (see [25, 26]).
Note that (1) degenerates when goes to zero. We transform the BlackScholes equations (1)(2) into a nondegenerate partial differential equation by using a log transformation where .
For applying the numerical method, we truncate the infinite domain into a truncated domain , where and are chosen properly so that, for practical purposes, they do not affect the option price. Based on Willmott et al.'s estimate [2] that the upper bound of the asset price is typically three or four times the strike price, it is reasonable for us to set , while, for the lower bound of the asset price, since is negative enough, we take for convenience in numerical experiments. Therefore, in the remaining of this paper, we will consider the following problem: Here, the right boundary condition is chosen according to Vázquez [8]. Normally, this truncation of the domain leads to a negligible error in the value of the option [25].
3. The Temporal Semidiscretization
To approximate the solution (4), first, we apply the implicit Euler method to discretize the temporal variable. This scheme, on a uniform mesh reads where and denotes the approximation of the exact solution at the time level .
Similarly to Kellogg and Tsan [27], we can prove that the differential operator satisfies a maximum principle, and, consequently, Hence, we can obtain the following result.
Lemma 1. The temporal semidiscretization scheme (6) is unconditionally stable.
Estimates for the global error are deduced from appropriate bounds of the local error, where the auxiliary problem is introduced to define the local error.
Lemma 2 (local error estimate). The local error associated with the method (9), defined by , satisfies
Proof. Using Taylor expansion, we have for . From (9) and (11), it is straightforward to show that the local error is the solution of the problem and, therefore, the result follows from the maximum principle for the operator .
Lemma 3 (global error estimate). The global error associated with the implicit Euler method (6), given by , satisfies and, therefore, the temporal semidiscretization scheme is a firstorder convergent scheme.
Proof. The global truncation error at time can be decomposed as By relations (6) and (9), we have Applying the maximum principle for the operator , we can obtain Thus, from (14)–(16), we have for , where is a positive constant independent from .
4. The Spatial Discretization
For the approximate solution of the semidiscretization problem (6), the spatial discretization is performed on a uniform mesh for the computational domain . Thus, at each time point , we apply a cubic spline scheme on the above uniform mesh to approximate problem (6).
Let be the approximate solution of the exact solution of the boundary value problem (6) at the th time level. At each subinterval , the cubic spline function has the following form: where , and are constants. Using the notation for approximation of at mesh points and as interpolatory constraints. From algebraic manipulation, we can obtain where . Using the continuity of the first derivative at mesh point , we get the following equation: Substituting in (21) and using the following approximations for firstorder derivative of we get the following spline difference scheme: where and .
It is easy to see that the matrix associated with the discrete operator is an Mmatrix for sufficiently small . Hence, the following discrete maximum principle holds true.
Lemma 4 (discrete maximum principle). For sufficiently small , the operator defined by (25) on the uniform mesh satisfies a discrete maximum principle; that is, if is a mesh function that satisfies , and , then , for all .
From the above lemma, we can conclude that the spatial discretization scheme (24) is maximumnorm stable.
To prove the convergence of the spline difference scheme, we discretize the auxiliary problem (9) and obtain
Lemma 5. Let be the solution of (9) and be the solution of (26). Then, we have the following error estimate:
Proof. We use a Taylor expansion at to obtain the following local truncation error estimate: Hence, using the discrete maximum principle (Lemma 4) for the discrete operator , we have which completes the proof.
5. The Fully Discrete Scheme
Combining the time semidiscretization scheme (6) with the spatial discretization scheme (24), we can obtain the following fully discretization scheme: where the discrete operators are described in Section 4 and is the fully discrete approximation to the exact solution of (4) at the mesh point .
Now, we can get the main result for our difference scheme.
Theorem 6. Let be the exact solution of (4) and let be the discrete solution of the fully discrete scheme (30). Then, the global error of our difference scheme satisfies where is a positive constant independent of and .
Proof. The global error at the time can be decomposed in the form From Lemmas 2 and 5, we can obtain Further, it is easy to see that can be written as the solution of one step of (30) with zero boundary conditions and as the final value. Applying the stability of the discrete scheme (Lemma 4), we have Then, from (33) and (34), a recurrence relation for the global errors follows, and, from it, the result of Theorem 6 can be obtained immediately.
6. Numerical Experiments
In this section, we present some numerical results to examine the performance and convergence of the cubic spline method. Errors and convergence rates for the numerical scheme are presented for three examples.
Example 1. European call option with parameters: , , , , , and : in this case, the analytical solution is
where
The maximum error is given in Table 1. The analytical and numerical solution profiles are given in Figure 1.

Example 2. European call option with parameters: , , , , and : as in Example 1, in this case, the analytical solution is known.
The maximum error is given in Table 1. The analytical and numerical solution profiles are given in Figure 2.
Example 3. European call option with parameters: , , , , , and : here, the volatility function is the same as the one given in Toivanen [28] and Kadalbajoo et al. [23, 24].
In this case, the exact solution is not known. We use the approximated solution of and as the exact solution. We present the error estimates for different and . Let denote “the exact solution." We measure the accuracy in the discrete maximum norm
and the convergence rate
The error estimates and convergence rates are listed in Table 3. The analytical and numerical solution profiles are given in Figure 3.
From Figures (1)–(3), it is seen that the numerical solutions by our method are nonoscillatory. From Tables 1, 2, and 3, we see that is close to 4, which supports the convergence estimate of Theorem 6. They indicate that the theoretical results are fairly sharp.


Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors would like to thank the anonymous referee for several suggestions for the improvement of this paper. The work was supported by the National Natural Science Foundation (Grant no. 11201430) of China, Ningbo Municipal Natural Science Foundation, Projects in Science and Technique of Ningbo Municipal (Grant no. 2012B82003) of China, and Key Research Center of Philosophy and Social Science of Zhejiang Province—Modern Port Service Industry and Creative Culture Research Center.
References
 F. Black and M. S. Scholes, “The pricing of options and corporate liabilities,” Journal of Political Economy, vol. 81, no. 3, pp. 637–654, 1973. View at: Google Scholar  Zentralblatt MATH
 P. Wilmott, J. Dewynne, and S. Howison, Option Pricing: Mathematical Models and Computation, Oxford Financial Press, Oxford, UK, 1993.
 J. C. Cox, S. Ross, and M. Rubinstein, “Option pricing: a simplified approach,” A Simplified Approach, vol. 7, no. 3, pp. 229–264, 1979. View at: Google Scholar  Zentralblatt MATH
 J. C. Hull and A. White, “The use of control variate technique in option pricing,” Journal of Financial and Quantitative Analysis, vol. 23, no. 3, pp. 237–251, 1988. View at: Google Scholar
 E. S. Schwartz, “The valuation of warrants: implementing a new approach,” Journal of Financial Economics, vol. 4, no. 1, pp. 79–93, 1977. View at: Google Scholar
 G. Courtadon, “A more accurate finite difference approximation for the valuation of options,” Journal of Financial and Quantitative Analysis, vol. 17, no. 5, pp. 697–703, 1982. View at: Google Scholar
 L. C. G. Rogers and D. Talay, Numercial Methods in Finance, Cambridge University Press, Cambridge, UK, 1997.
 C. Vázquez, “An upwind numerical approach for an American and European option pricing model,” Applied Mathematics and Computation, vol. 97, no. 23, pp. 273–286, 1998. View at: Publisher Site  Google Scholar  MathSciNet
 Z. Cen and A. Le, “A robust finite difference scheme for pricing American put options with singularityseparating method,” Numerical Algorithms, vol. 53, no. 4, pp. 497–510, 2010. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 Z. Cen and A. Le, “A robust and accurate finite difference method for a generalized BlackScholes equation,” Journal of Computational and Applied Mathematics, vol. 235, no. 13, pp. 3728–2733, 2011. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 Z. Cen, A. Le, and A. Xu, “A secondorder difference scheme for the penalized blackscholes equation governing American put option pricing,” Computational Economics, vol. 40, no. 1, pp. 49–62, 2012. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 S. Wang, “A novel fitted finite volume method for the BlackScholes equation governing option pricing,” IMA Journal of Numerical Analysis, vol. 24, no. 4, pp. 699–720, 2004. View at: Publisher Site  Google Scholar  MathSciNet
 L. Angermann and S. Wang, “Convergence of a fitted finite volume method for the penalized BlackScholes equation governing European and American option pricing,” Numerische Mathematik, vol. 106, no. 1, pp. 1–40, 2007. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 N. Rambeerich, D. Y. Tangman, A. Gopaul, and M. Bhuruth, “Exponential time integration for fast finite element solutions of some financial engineering problems,” Journal of Computational and Applied Mathematics, vol. 224, no. 2, pp. 668–678, 2009. View at: Publisher Site  Google Scholar  MathSciNet
 K. M. Liew, Z. X. Lei, J. L. Yu, and L. W. Zhang, “Postbuckling of carbon nanotubereinforced functionally graded cylindrical panels under axial compression using a meshless approach,” Computer Methods in Applied Mechanics and Engineering, vol. 268, pp. 1–17, 2014. View at: Publisher Site  Google Scholar  MathSciNet
 L. W. Zhang, P. Zhu, and K. M. Liew, “Thermal buckling of functionally graded plates using a local Kriging meshless method,” Composite Structures, vol. 108, pp. 472–492, 2014. View at: Google Scholar
 R. J. Cheng, L. W. Zhang, and K. M. Liew, “Modeling of biological population problems using the elementfree kpRitz method,” Applied Mathematics and Computation, vol. 227, pp. 274–290, 2014. View at: Publisher Site  Google Scholar  MathSciNet
 L. W. Zhang, Z. X. Lei, K. M. Liew, and J. L. Yu, “Static and dynamic of carbon nanotube reinforced functionally graded cylindrical panels,” Composite Structures, vol. 111, pp. 205–212, 2014. View at: Google Scholar
 L. W. Zhang, Y. J. Deng, and K. M. Liew, “An improved elementfree Galerkin method for numerical modeling of the biological population problems,” Engineering Analysis with Boundary Elements, vol. 40, pp. 181–188, 2014. View at: Publisher Site  Google Scholar  MathSciNet
 C. C. Christara, T. Chen, and D. M. Dang, “Quadratic spline collocation for onedimensional linear parabolic partial differential equations,” Numerical Algorithms, vol. 53, no. 4, pp. 511–553, 2010. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 M. Holtz and A. Kunoth, “Bsplinebased monotone multigrid methods,” SIAM Journal on Numerical Analysis, vol. 45, no. 3, pp. 1175–1199, 2007. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 M. H. M. Khabir and K. C. Patidar, “Spline approximation method to solve an option pricing problem,” Journal of Difference Equations and Applications, vol. 18, no. 11, pp. 1801–1816, 2012. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 M. K. Kadalbajoo and L. P. Tripathi, “A robust nonuniform Bspline collocation method for solving the generalized BlackScholes equation,” IMA Journal of Numerical Analysis, vol. 34, no. 1, pp. 252–278, 2014. View at: Publisher Site  Google Scholar
 M. K. Kadalbajoo, L. P. Tripathi, and A. Kumar, “A cubic Bspline collocation method for a numerical solution of the generalized BlackScholes equation,” Mathematical and Computer Modelling, vol. 55, no. 34, pp. 1483–1505, 2012. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 R. Kangro and R. Nicolaides, “Far field boundary conditions for BlackScholes equations,” SIAM Journal on Numerical Analysis, vol. 38, no. 4, pp. 1357–1368, 2000. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 O. A. Ladyzenskaja, V. A. Solonnikov, and N. N. Ural'ceva, Linear and Quasilinear Equations of Parabolic Type, American Mathematical Society, Providence, RI, USA, 1968. View at: MathSciNet
 R. B. Kellogg and A. Tsan, “Analysis of some difference approximations for a singular perturbation problem without turning points,” Mathematics of Computation, vol. 32, no. 144, pp. 1025–1039, 1978. View at: Google Scholar  Zentralblatt MATH  MathSciNet
 J. Toivanen, “Numerical valuation of European and American options under Kou's jumpdiffusion model,” SIAM Journal on Scientific Computing, vol. 30, no. 4, pp. 1949–1970, 2008. View at: Publisher Site  Google Scholar  MathSciNet
Copyright
Copyright © 2014 Jian Huang and Zhongdi Cen. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.