Introducing the Power Series Method to Numerically Approximate Contingent Claim Partial Differential Equations

We introduce a previously unused numerical framework for estimating the Black-Scholes partial differential equation. The approach, known as the Power Series Method (PSM), offers several advantages over traditional finite difference methods. Our objective is to highlight the advantages of the PSM over traditionally used numerical approximation approaches. To meet this we deploy a numerical approximation scheme to illustrate the PSM. The PSM is more stable than explicit methods and thus computationally more efficient. It is as accurate as hybrid approaches like Crank Nicolson and faster to compute. It is more accurate over a far wider spectrum of time steps. Finally, and importantly, it can be expressed analytically thus offering the capability of performing comparative statics in a far more stable and accurate environment. For a more complex application this last advantage may have wide im-plications in producing hedge ratios for synthetic replication purposes.


Introduction
Numerical approaches to solving contingent claim based partial differential equations (PDE) have focused on finite difference method (FDM) approaches. Specifically, these have centered on explicit FDM (EFDM) approaches and Crank-Nicholson methodologies (CNM) and modification of these methods. Here we consider the approach introduced by Parker and Sochacki (1996,2000) [1] [2] that uses a power series approximation based upon the Picard method.
The former methods provide solutions at grid points determined by the user. The method introduced by Parker and Sochacki provides a function that is defined over an entire time interval and is an accurate approximation to the solution to the contingent claim problem. Since one has a function that approximates the solution accurately, one can differentiate, integrate or use this function to solve at what time a certain value is achieved. These advantages become more apparent as the contingent claim becomes more complicated. The objective here is to illustrate the advantages relative to plain vanilla options in an effort to present why this approach may be a superior framework for complicated, nonlinear, or less stable contingent claim PDE's. In the literature this method is also known as automatic differentiation (Gofen [3] (2009) and Neidinger [4] (2010)), differential transform method (Mirzaee [5] (2011)), the Parker-Sochacki method (Stewart and Bair [6] (2009) and Rudmin [7] (1998)) and the power series method. In this article, we will use the terminology power series method (PSM) for the method that was discovered by Parker and Sochacki through Picard iteration.
Our paper follows Anwar and Andallah [8] (2018) and Cen and Le [9] (2011) by introducing PSM as an alternative and improved numerical scheme to estimate the Black Scholes PDE. The PSM approach offers advantages over traditional FDM that produces both operational and estimation improvements. Moreover, it enables the user to better approximate solutions to the PDE in far greater levels of precision. Also, PSM does not have to be limited to an FD mesh. We demonstrate this below in some examples that highlight these points. The highlights will be emphasized in our examples. We also point out that PSM is easy to program in both Matlab and Maple and is a natural extension of EFDM. In future papers we will show that it works for variable rates, dividends and volatility including those that produce a nonlinear PDE.
Since this method is not common in the mathematical financial community, we will demonstrate it through some examples and compare it to EFDM and CNM. As the reader goes through the examples, some important concepts should become clearer. One is that we are determining the coefficients of the power series solution to the differential equations by two approaches. One method is simpler and easier to program than the other. However, the second method is also valid and gives us the ability to do mathematical analysis on the differential equation, the solution and the properties of the approximate solution. By having a polynomial estimate, we are able to analyze the discontinuity in the derivative at the strike price and study oscillations in difference methods noted by many researchers, including Cen and Le [9]. One can also generate comparative statics commonly referred to as the Greeks using calculus. The polynomial is easily differentiated at any combination of underlying variables and thus offers a framework where analytics on the contingent claim valuation is greatly facilitated.
In this article, our goal is to show that the method is a powerful method for obtaining approximate solutions to contingent claims and is superior to discrete

Introduction of the PSM Approach to Estimating Differential Equations
In this section, we introduce the general approach of PSM to differential equations. This is necessary in order to develop the approach as it specifically applies to the contingent claim partial differential equation presented in section 3 of the paper. To introduce the framework we start with a simple asset, V, that increases over time (t) due to continuous compounding at a constant rate (r). The ordinary differential equation (ODE) describing this process is where V 0 is the initial value (IV) or initial condition. The solution to this ODE is Note that the solution depends on V 0 , r and t. That is, we can write to highlight the dependence of the value of the solution on V 0 , r and t.
In EFDM and CNM one discretizes time to approximate the solution to this problem. That is, one chooses specific time values at which an approximation for V will be obtained. We assume these times are n t n t = ∆ for 1, 2,3, , The CNM approximation for IV ODE (1) is obtained from ( )  It is also well known that by choosing a large enough counting number J that 0 J j j j V t = ∑ will be as accurate an approximation as one desires to That is, one has a function of that accurately approximates the true solution We now demonstrate two methods for obtaining the PS form for The second method is to substitute Equating like powers of t in this last equation gives the iterative sequence Writing out a few of the iterates gives 1 and so on. Again, one can determine that these are the coefficients of the PS for the solution to the ODE (1). The recurrence relation (2) is easily programmed in either a symbolic or numerical computing environment. Therefore, one can, in principle, generate as accurate a polynomial solution to the true solution as desired by generating enough coefficients using recurrence relation (2). That is, one code, as will be demonstrated in this article, can generate as high an order of accuracy as one desires.
One important item to note is that if 1 J = then the approximation for ( ) V t using PSM is given by V rV t + ∆ which we noted above is the EFDM approximation. This is an important result as we now can generate a single framework that includes PSM and EFDM as a special case ( 1 J = ). If 2 J = , one obtains as the approximation of ( ) V t over the interval 0 t t ≤ ≤ ∆ . To extend to the time interval over the interval 2 t t t ∆ ≤ ≤ ∆ , one updates 0 V to be ( ) and then generates a new 1 2 , V V with this 0 V and obtains an approximation for ( ) . One continues for any desired N. Of course, the process is similar for any J. The larger J is or the smaller t ∆ is the better the approximation will be. Warne P.G. et al. [10] use the Picard iterates and PS to derive an a-priori error bound for PSM using this fact. This is how a PSM approximation for a plain vanilla option will be derived in our numerical analysis below.
It is important to mention a few other points. Picard showed that his iterative method works for a large class of ODEs. That is, the Picard iterates will converge to the solution to the ODE for many ODEs. This is important for our purposes and why offering the Picard approach as a polynomial solution. Parker and Sochacki [1] used this fact to generate solutions to nonlinear ODEs. In fact, the second method also works for nonlinear ODEs which was shown by Parker and Sochacki and exploited to give convergence and error results by Carothers et al. [10] [11] [12]. The important point is that we now have two methods for developing algorithms and proving theorems for solutions to linear and nonlinear ODEs. Also we point out that since the two methods are computational and based off the ODE, they can compute either a symbolic solution which can be analyzed in terms of the parameters defining the ODE, like V 0 and r in ODE (1) or strictly numerically for visualization and numerical analysis. This offers tremendous advantages over the traditional FDM approaches as it enables us to work with the actual solution of the PDE not only discrete numerical approximations of the PDE. This is not a subtle distinction but a meaningful shift in how to understand the valuation dynamics of contingent claims. Now we can actually operate on a time continuum and almost any space dimension required with more accuracy. By operating in this framework we are able to more accurately understand the PDE dynamics and corresponding first and second order hedging or replication applications. This will be demonstrated for plain vanilla options using the Black-Scholes model in this article and on various other options, including nonlinear options in future papers. We now demonstrate the two methods on a nonlinear ODE in order to indicate how we will use PSM on nonlinear options in our future papers.
We modify the IV ODE (1) to for r, l positive. This dynamical system puts a bound on the growth of ( ) Determining the PS form of ( ) 0 ; , , V t r l V is not an easy task. However, if one uses the two methods presented in our first example, it is straightforward. First, we note that if we have two PS then the PS product of A and B is given by Cauchy's formula We note that the Picard iterates are much more complicated than the Picard iterates in the linear IV ODE of our first example, but that a pattern is arising. Even for this complicated example, Picard showed that Therefore, using Cauchy's formula we have that ( ) This gives the recurrence relation We leave it to the reader to verify that 1 2 3 , , V V V given by this recurrence relation agrees with the coefficients in the Picard iterate [ ] 3 P above. Therefore, the solution to the IV ODE (3) is given by Again, we can approximate ( ) V t to any accuracy desired over an interval 0 t t ≤ ≤ ∆ for some chosen t ∆ and counting number J using on that interval. Note that the coefficients, once again, depend on V 0 , r and l.

Black Scholes Plain Vanilla Option
Using the development in section 2 we can now turn our attention to option valuation. In a plain vanilla option, the value of the option, V depends on the price of an underlying asset, S and time, t. Therefore, one wants to determine where V represents the value of the option, σ the standard deviation of returns on the underlying asset, S, and r is the risk-free borrowing and lending costs over the life of the contingent claim, and q is cash flow emanating from S. All rates are continuous. F is the value of the PDE at any given combination of asset value or time and is zero for simple contingent claims. In this article, we will assume F is zero, but in a future paper we will consider arbitrary F.
We will compare the EFDM, CNM and PSM for IV PDE (4) to the closed form solution. Again, since PSM gives us a function of t, we can compare PSM to the closed form solution at any t for a given i S .
As for the two IV ODEs above, Equation (4) can be written equivalently as One can then build a Picard iterative process for V using this IE.  , S τ of interest. Now that we have this solution, we can analyze the dependence of V on not only on S, but also σ, q, r and K.
Note that this method will generate a power series solution (polynomial) approximation for V for more complicated situations, including a polynomial F, an r that depends on S and/or t or even V and other Q. Thus, this approach can be used for the aforementioned nonlinearities. This allows one to do a thorough parameter analysis on K, r, q, σ, F and Q in a very wide-ranging set of applications. We will leave the more complex applications to future research. Here, we want to introduce the approach and compare it to the aforementioned common numerical approaches found in the literature.
Since the EFDM and CNM are discrete methods, we let and we let n i W be the approximation to ( ) i v n t ∆ . We mention that PSM can be used in the non-discretized case as shown above and an example generated in Maple will be given below. The discretized version of Equation (4) that we consider is then Equation (7) in discretized form becomes the ODE with the initial conditions for some counting number I, we have a system of I IV ODEs that we can solve using PSM.
That is, we assume for each I and some user chosen J. The larger J is the more accurate an approximation ( ) i v τ will be to ( ) , i V S τ . We can use the PSM in either a symbolic or numeric computing environment. This bypasses the restriction of using only time discrete approaches like EFDM and CNM. For our purposes here, we will contrast the PSM approach to these two numerical approaches to illustrate its strengths. We will use only the PS approach for PSM because of its ease to program. Since the first PSM polynomial (PSM 1) was shown to be equivalent to EFDM, we can generate EFDM with our PSM code. The Picard method is a useful tool for analysis, The PSM effectively offers a more robust framework for contingent claim valuation and a corresponding analysis. 1 We will demonstrate that the PSM framework is a powerful setting to be applied to numerically understanding contingent claim valuations and more im- 1 In other words, Equation (6) can be used to create polynomial estimates to the solution to Equation (4) over the entire valuation spectrum. These polynomials are then easily manipulated to produce estimates to all of the Greeks throughout the same spectrum. These will be more accurate than any FDM approach and also require a single computation. FDM requires computations over several input variations to obtain the Greeks or any kind of comparative static analysis. These repeated numerical computations may introduce a propagation of errors and may result in a compounding error. PSM does not suffer from this issue.
portantly the corresponding comparative statics relationships often used in practice by risk managers and traders. This is an important development since most contingent claim-based valuation models do not have analytical solutions so having another numerical approach like PSM may prove invaluable.
The EFDM approximation to Equation (8) that will be used is ( ) This form shows the slope term multiplied by τ ∆ .
The CNM approximation to Equation (8) We solve this system of linear equations using a tri-diagonal matrix solver.
We compare EFDM, CNM and PSM for polynomials of various degrees with the closed form solution to IV PDE (4). For PSM we will determine an accurate PS solution for the ODEs in Equation (8). We will demonstrate the power of being able to choose J with- Again, we point out that one does not have to discretize in S to use PSM as was shown in our third example. We are doing this in this article for comparison with EFDM and CNM. We will also present some of the coefficients that PSM produces so that one can view the approximate polynomial solutions to the "true" answer. We also point out that polynomials are straightforward to evaluate efficiently on a computer using Horner's algorithm which we actually do here.

Motivation behind PSM
To our knowledge the PSM approach has not been extensively used within the finance literature, if at all. The advantages of PSM are many and several of these are particularly important in the contingent claim valuation framework. The neuroscientists Stewart and Bair [6] showed that PSM was competitive with the Runge-Kutta order 4 method (RK4) and the Bulirsch-Stoer method (BSM) and, in most cases, superior for Hodgkin-Huxley type differential equations, which are nonlinear. The astrophysicists Pruett, Rudmin and Lacy [35] showed that PSM was in many cases superior to RK4 and BSM for Newton's N-Body problem. Also, the astrophysics teams Nurminskii and Buryi [36] and Pruett, Ingham and Herman [37] and the neuroscientists Synkiewicz [38] and Yudanov, Shaaban, Melton and Reznik [39] have shown that it is in general easier to parallelize PSM codes and obtain close to linear speed up than other codes, including on graphical processing units.
As outlined above the PSM can be used symbolically and/or numerically to analyze or obtain numerical results. A symbolic polynomial estimate to the PDE solution directly offers the user the ability to compute not only the value of the contingent claim but also the corresponding comparative statics of that value to changes in the underlying inputs. These are well known as the Greeks in the finance literature. Once the symbolic solution is created no further numerical estimation is required within PSM. This is very different that any FDM approach.
All FDM approaches would require repeated numerical estimations to compute the Greeks. The PSM doesn't require that and so errors introduced are not propagated further like in the FDM framework. We will illustrate this below.
The numerical PSM approach also subsumes the EFDM as a special case and thus can always be at least as accurate and stable as that approach. This enables the user to create a single PSM program and be able to easily use it when the EFDM is needed. Additionally, for similar reasons, the stability requirement of the PSM is never more constraining than the EFDM approach.
The CNM approach has the advantage of most always being stable but it is particularly inefficient computationally. CNM will also not have solutions for more complicated valuation applications that result in nonlinear frameworks.
PSM is both more computationally efficient and in almost all cases more accurate. It also lends itself to far more complicated problems than the CNM framework.
Computational errors within the PSM framework also tend to remain of the same order regardless of the underlying values of the inputs whereas those associated with FDM can compound (grow) throughout the FD grid.
Money [40] showed that when applied discretely to PDEs, PSM is a generalization of the Lax-Wendroff method (LWM) that is commonly used and/or modified in computational fluid dynamics that involve nonlinearities. He determined stability conditions for some linear PDEs, including the heat equation. His work was for constant coefficient PDEs, but he did show its applicability to some nonlinear PDES. As far as we know, this is the first application of PSM to PDEs with non-constant coefficients. In the Black-Scholes-Merton options modeling PDE (BSMOMPDE) the coefficients depend on the variable being discretized. In future papers, we will demonstrate the relationship between PSM and LWM in both the linear and nonlinear case and derive stability and convergence conditions. Here our goal is to introduce the method, to show its efficacy and its potential strengths for variable coefficient PDEs and show that we can improve on the stability of the EFDM and be competitive with CNM. It is not an easy matter to extend CNM to nonlinear BSMOMPDEs. However, PSM extends naturally without having to solve systems of nonlinear equations to get a numerical answer.
For the discretizations we used above the EFDM for the BSMOPMPDE has an error that is The CNM for the BSMOPMPDE has an error that is which is an improvement over the explicit FDM.
The PSM of order k (using polynomials of degree k) has an error that is which allows analyzing the error in using PSM on the BSMOPMPDE.
Finally, even the analytic solutions to the contingent claim BSMOPMPDE requires numerical estimation since it contains the cumulative distribution function. The distribution function requires a numerical approximation to the integration. This also introduces an error. It is possible that the PSM may more accurately approximate the true solution than this approximation of the analytic solution. Also, since one can increase the order of PSM, one has the advantage of studying convergence with respect to the grid size, the order of the polynomial used in PSM and the time step. This allows the potential of the computer giving as accurate an answer as possible.   This is the rare case where EFDM is the most accurate. Somewhat counter intuitively the increased terms in the PSM introduce a source of small round off errors when the time step is very small. However, since the EFDM is just a special case of the PSM we can determine easily at what time step does the PSM overtake the EFDM in accuracy. That is an important determination because as the time step increases the speed and stability of the PSM will dominate the tradeoff quite quickly. To illustrate this we increase the time step in the next example. Here we can plainly see that the EFDM is far less accurate than all the other methods and begins to oscillate. This is due to the instability of the EFDM. The second order PSM (PSM 2) and fourth order (PSM 4) is as accurate as the CNM and everywhere more accurate than the EFDM. As the time step increases and EFDM becomes less stable, and loses accuracy. PSM 4 is as accurate as CNM and close to expiration sometimes more accurate. This implies that the PSM can be computed more quickly than both CNM and EFDM without losing accuracy. This is particularly important in a trading environment or any application where resources are constrained. As we extend our research to non-linear applications this property will become even more important.

A Note on Discontinuity
Cen and Le (2011) highlight many of the difficulties within numerical routines in dealing with the non-differentiability at the strike price. In Figure 2 we can see the oscillation aspects to some degree within the EFDM framework. The general symbolic PSM method cannot be used on a non-differentiable initial condition. Here we approximate the put payoff with an infinitely differentiable function which gives us advantages over Cen and Le [9]. This is a relatively common issue throughout continuous mathematics. Whenever an application contains a discrete boundary condition difficulties will arise by the very nature of differentiability. Consequently, numerical approaches like FDM and PSM will have to deal with this issue. As we illustrate below the PSM easily addresses this problem that can negatively impact the FDM approaches.
One reason that PSM does not work in a symbolic environment for BSMOPMPDE is that the payoff functions are not differentiable at the strike price, K since the payoff for a call is given by ( ) ( ) max , 0 Q S S K = − and for a put is given by ( ) ( ) max , 0 Q S K S = − and to do PSM on a PDE the initial condition must be infinitely differentiable. In the Appendix, we give the second degree polynomial approximation to ( ) , V S τ that Maple generated (The fourth degree polynomial is too large to represent.) and the infinitely differentiable Q we used to approximately the put payoff even at the discrete boundary condition.

Conclusion
We have introduced the PSM framework to the finance literature and illustrated its advantages over traditional FDM. It is more accurate and offers a polynomial representation which can be further leveraged for both higher accuracy and more computational efficiencies. These advantages were illustrated through numerical examples. Specifically, as the time step increases (and computational speed increases) the advantages of PSM over EFDM becomes more apparent. Additionally, the PSM is everywhere as accurate as the CNM without the computational limitations. The PSM easily addresses discrete boundary conditions as well. Finally, due to the polynomial representation of PSM, analysts can far more accurately approximate values and hedging parameters for almost any set of inputs without being constrained by a finite difference grid. Future research will include more complicated nonlinear applications that will highlight the advantages over traditional FDM even more. These future applications offer fertile ground for research. Both FDM and CNM approaches struggle with nonlinear applications. PSM offers a framework that doesn't suffer the same limitations. PSM is an excellent alternative to the numerical finance literature.