Using the Power Series Method to Evaluate Non-Linear Contingent Claim Partial Differential Equations ()
1. Introduction
Buetow and Sochacki [1] use the PSM approach to evaluate linear PDE generated within the contingent claim valuation framework. Here, we extend their approach to nonlinear PDE (NLPDE). NLPDE occur in contingent claim valuation schemes under various settings. Leland [2] introduced transactional inefficiencies to the original Merton [3] and Black and Scholes [4] frameworks which resulted in a NLPDE. Several variations followed as did approaches to try and solve them. Traditional numerical schemes that solve contingent claim partial differential equations (PDEs) abound in the literature. Our focus here is to extend Buetow and Sochacki [1] into the NLPDE framework; their paper is the first to introduce the Parker and Sochacki [5] [6] Power Series Method (PSM) method within the contingent claim literature.
A contingent claim is a financial asset that has a future payoff contingent on an underlying asset meeting certain criterion. In this paper, we focus on options that have a payoff contingent on the position of the underlying asset relative to a strike price. Traditional valuation frameworks assume that the underlying asset is infinitely liquid; that is, it can be traded with no transaction costs, in infinitesimally small increments, and on a continuum. This assumption is clearly not realistic. Removing it from the valuation process results in altering the valuation framework from a linear PDE to a NLPDE. We will focus on the impact of the liquidity parameter on the valuation estimates to highlight the benefits of using PSM over alternative numerical approaches.
Buetow and Sochacki [1] presented the advantages of the PSM within the contingent claim valuation framework by approximating the plain vanilla PDEs. They illustrate that the PSM approach is inherently more stable than traditional finite difference approaches and so more computationally efficient. PSM is also as accurate as hybrid finite difference approaches like Crank Nicolson and computes a numerical solution more quickly. PSM also is more accurate over a far wider spectrum of time steps than alternative numerical approximation frameworks. Additionally, PSM can be expressed analytically thus enabling the computation of comparative statics (i.e., Greeks) within a far more stable and more accurate environment. Here, we extend their analysis to NLPDEs.
Examples of NLPDE’s in the literature include Amster and Mogni [7], Amster, Averbuj, Mariani, and Rial [8], Bakstein, Howison [9], Barles and Soner [10], Boyle and Vorst [11], Ehrhardt [12], Frey [13], Frey and Patie [14], Frey and Stremme [15], Wilmott [16], and Hoggard, Whaley, and Wilmott [17]. While this is not an exhaustive list, it certainly represents a wide spectrum of frameworks within which a NLPDE results to evaluate contingent claims under various environments where the assumptions of the original Black and Scholes and Merton are relaxed. For our development here we will focus on the NLPDE resulting from transactional inefficiencies as developed by Schonbucher and Wilmott [18], Frey [19], and Liu and Yong [20]. Transactional inefficiencies like illiquidity and trading costs are certainly the most practical. All market participants, both large and small, face these realities when trading, hedging, or replicating contingent claims. Hence, our focus is on this NLPDE.
The aforementioned references focused on developing the necessary mathematics ensuring that the framework met the necessary principles for a solution. Many of these require numerical approaches to enable a solution to be quantified. Other numerical and analytic solutions can be found in Abe and Ishimura [21], Ankudinova and Ehrhardt [22], Anwar and Andallah [23], Arenas, Gonzalez-Parra, and Caraballo [24], Cen and Le [25], Company, Jodar, and Pintos [26], Cont and Voltchkova [27], Duffie [28], Forsyth and Labahn [29], Giles and Carter [30], Gonzalez-Parra, Arenas, and Chen-Charpentier [31], Jandacka and Sevcovic [32], Koleva and Vulkov [33], Kutik and Mikula [34], Lesmana and Wang [20] [35], Sevcovic [36], Sevcovic, Stehlıkova, and Mikula [37], Sevcovic and Zitnanska [38], Wang and Forsyth [39], Zhou, Han, Li, Zhang, and Han [40], Zhou, Li, Wei, and Wen [41], Brennan and Schwartz [42], and Buetow and Sochacki [43] [44] [45]. A good analysis of the Frey model can be found in [19] and Polte [46]. The methodologies and analysis in the above studies often require various transformations into linear approximations of the original NLPDE to be implemented. Our approach here offers what we believe to be a more robust framework for evaluating the NLPDE in the transactional inefficiency environment.
2. Model
Buetow and Sochacki [1] demonstrated that PSM can be used to efficiently solve and analyze the forward linear BS options model and is comparable to the explicit finite difference method and the Crank-Nicolson method. Matić, Radoičić and Stefanica [47] show that PSM can be used to efficiently solve the inverse problem for determining the implied volatility in the linear BS options model. In this paper, we will present how to use PSM to solve and analyze the forward nonlinear BS options model. To be specific, we choose Frey’s model. However, the method presented can easily be extended to other nonlinear BS models.
As in [47], we use auxiliary variables and power series to solve this nonlinear BS options model. In the appendix, we give the symbolic solution for this model similarly to how it was done in [1]. As in [1] and Glover, Duck, Newton [48], we approximate the pay off with a smooth function. We do it to generate a symbolic solution that can be analyzed using calculus including for the “Greeks”. However, since we can also generate polynomials that approximate the solution, one can also use these polynomials to do the same analysis. The symbolic polynomials presented in the appendix were done in the software package Maple.
In [1], we derived methods to approximate the solution to the linear BSOPM using finite differences incorporated with the method of lines to develop a system of ODEs. We then used PSM to approximate the solution to this system of ODEs. This is basically an extension of the Lax-Wendroff method for numerically solving PDEs to higher derivatives. With PSM, this is done automatically (and symbolically) through the computer. In that paper, it was also demonstrated that if the payoff function is approximated with a smooth function one can generate a symbolic solution to the linear BSOPM PDE. In this paper, it will be demonstrated that this can be done for the nonlinear BSOPM PDE through the nonlinearity introduced by [19] and others and developed further by Wilmott. It is very useful that one can generate a symbolic solution to this NLPDE for analysis and numerical evaluation. Evidence of this is presented in the appendix using Maple.
3. Mathematical Framework
Relying on the seminal work of Kyle [49] and employing the classic no arbitrage replication framework commonly found within the option pricing literature, the following NLPDE is easily derived:
(1)
where the diffusion term U implemented in this paper, which represents the transactional inefficiencies will be
This term introduces the nonlinearity. The value of the option is given by
where S is the value of the underlying asset for the option and t is calendar time. We denote our payoff as
where T is the maturity time of the option.
Equation (1) has appeared extensively in the literature. Most commonly, the exponent, a, in
is set to −2. We will use that quantification in this work. [18] assume that the function
is a dimensionless constant measuring the liquidity of the market. [19] has a similar form for
. [20] define
as
where
is a coefficient and a constant that measures the price impact,
is time to expiry and
is a decay coefficient. We follow [15] [18] and [19] and assume that
is a constant. [48] present Equation (1) with constant
but include additional insights into its various properties. There is wide acceptance of this form with the most practical implications throughout the industry.
In the linear paper ( [1]), the put1 payoff function
(where Kis the strike price of the option) was approximated by the smooth function
. In [48], the put payoff is approximated by
. Letting
gives a better approximation. Unfortunately, this function is not smooth at
. Therefore, we approximate the put payoff here by
which leads to a better approximation of the put payoff as
. All the numerical examples presented in this paper were done using
and
with
and gave an answer to within 10−5 of the results presented.
1For a call,
is replaced by
and
is replaced by
in each
for
.
As shown in the linear paper, one can obtain a symbolic answer using a smooth Q. There we showed how Picard iteration on the integral form of the PDE can generate the power series solution in t or
for the option value V. The same can be done for any of the nonlinear forms for the BSOPM. We provide the PSM outline for doing this for the nonlinear problem done in this work in the appendix where a symbolic solution up to degree 2 in
for any smooth Q, including
is provided. One can use this symbolic form of the solution to analyze the “greeks” of the option. We leave this as an exercise for the reader. In this paper, only the discretized PSM presented in the linear paper will be applied to the nonlinear BSOPM. In the appendix, we give a system of PDEs that can be used to obtain a symbolic (or numeric) power series solution to
for any smooth Q.
We also point out that since we have a power series in
, one can replace
by
(or
by a power series
) for any function z to get an expansion of V in powers of z. This allows one to get numerical approximations of V in terms of z or symbolic approximations of V in terms of z.
We will do our numerical derivation in
and consider only the following PDE
(2)
where
replaces
in Equation (1). We note again that our derivation can be used for more than the
given above.
4. Setting up the DPSM
We now outline how to obtain a Lax-Wendroff approximation to any degree desired in
for Equation (2) with the
given above. Hoffman and Frankel [50] and Evans, et al. [51] are both good and thorough introductions to Crank-Nicolson, MOL and Lax Wendroff schemes for linear and nonlinear PDEs. Money [52] presents a thorough demonstration of how PSM with MOL gives a generalization of the Lax Wendroff scheme with error estimates and stability analysis.
We will assume
is a constant, but the discretized PSM applies for most forms for
in the literature and certainly those discussed above. The derivation would be similar as presented below but with different auxiliary variables. A thorough discussion of auxiliary variables is given in Parker and Sochacki [5] [6] and in Sochacki [53].
We first deal with
. Note that if x is a power series in z and
then differentiating leads to
and then
(3)
Using this, it is straightforward to develop the power series coefficients for y. In fact, if
then
and substituting these four power series into Equation (3) gives
Multiplying out the series using Cauchy products as outlined in the linear paper produces
Equating the terms multiplied by like powers of z, and solving for the coefficient
in terms of all the lower coefficients of y provides the recurrence relation for each coefficient of the power series for y. Doing this, gives us
and then
Rewriting slightly gives us
Solving on the left-hand side for
leads to
or after re-indexing
(4)
This recurrence relation gives all the coefficients of the power series for
. Since
and PSM gives the power series for
in
, we can use the above recurrence relation to obtain the power series for U for any a but in our numerical examples we use
. We choose
because this is the result that occurs in the derivation using Ito’s lemma as shown by [19] and many others.
If
then
.
Since
for
, one has that
for
. If
is substituted into Equation (2), the PDE in Equation (2) gives rise to
One can rearrange this last PDE to read as
where the nonlinear terms can now be viewed as
as in Equation (1). This equation shows how sensitive the nonlinearity can become as x gets further from 0 (but stays < 1 in magnitude). (If
, one obtains the linear Black Scholes contingent claim PDE.) This sensitive nature created by the Frey illiquidity has been noted by other researchers (Heider [54]). However, we believe this is the first paper that studies this nonlinearity in a power series expansion as presented in this last equality. Some of the findings are not intuitive and show the complicated convergence (and stability) regions that arise in numerically solving this nonlinear PDE. We emphasize that the method presented works for
. However, our study shows that the size of
certainly determines the error, convergence, and stability in the numerical approximations. It probably also reaches to the limits of the nonlinearity to model a real economy.
5. The Discretization Scheme Using PSM
We now layout the finite difference grid for
. This is the same grid as given in the linear paper. We present it here for completeness.
We choose an increment value ∆S for the interval
and let
for
so that
and
. We determine a power series for
for each i. To do this, we let
be our Maclaurin polynomial approximation to
and let
be the value
. The discretized version of Equation (2) is
(5)
with
,
, D a discretized approximation for
and
a discretized approximation for
. For example, if we use the centered difference approximations
then Equation (5) in discretized form becomes the ODE
(6)
with
.
Therefore, if
for some counting number I, we have a system of I IV ODEs that we can solve using PSM. That is, we assume
(7)
for each i and some user chosen J that determines the degree of the Maclaurin polynomial that approximates
. We determine the polynomial for
using the recurrence relation (4) and the fact that we know the polynomials for
. The standard conditions for numerical solutions to Equation (2) (these are the one used in the linear paper) are used to give values to
and
. Since the PSM polynomial of degree
(PSM 1) is equivalent to EFDM, this PSM algorithm can be used to compare EFDM against PSM J for any user chosen J. It is straightforward to combine the right-hand side of Equation (6) into a power series in
after using Cauchy product of two power series to multiply out
.
We can compute the
numerically or symbolically using a software package. We present symbolic representations using Maple in the Appendix. Maple was able to generate an 8th degree symbolic polynomial for each
in less than two minutes. It did this also for
but the output is too large to present in this paper. In the Appendix the output is shown for the solution to Equation (2) with
using Picard iteration in Maple. It is also important to point out that one can replace as many as parameters with numerical values as desired in this symbolic polynomial to study the impact of other parameters in the Greeks. We note that coefficient of
has the term
in it and that this term occurs to more negative exponents in the coefficients of
with larger exponents along with terms in
. These contribute to the nonlinear feedback in an illiquid market.
Of course, one can use Equation (7) and discrete approximations for S derivatives to get the greeks numerically to order J.
Using
and
the fourth degree (Maclaurin polynomial) approximation to
(
) for a call at
is
.
Again, it is straightforward for Maple or a computing language to generate higher degree polynomial approximations with the same code. The only number the user must change is J.
In our numerical results, we demonstrate that the explicit FDM gives similar results to PSM under convergence regions but that PSM for
has larger (stronger) convergence regions. Most of the results shown with PSM are with
but we have verified these results with
. Since PSM can record any (all) of the Greeks, one can measure any of these quantities during the computing of the solution to determine their impact (size) on the value of the option and/or the error, convergence, stability of the numerical approximations.
6. Modeling Results and Comparisons
2Again, one should note that this payoff function is not everywhere differentiable either and so better alternatives exist. We offer such an approach in the appendix. For our purposes here we are focusing our presentation on a numerical approach such that the concerns associated with the payoff functional form are mitigated within the framework.
Figure 1 and Figure 2 clearly show the instability issues using the traditional explicit FDM approaches as also presented by [48]. The magnitude of the pricing impacts is extraordinarily large and oscillatory. The efficacy of the approach is clearly challenged. Figure 1 illustrates pricing impact oscillations an order of magnitude that are 104 the size of the option prices. Figure 2 shows that the estimated pricing impacts from the feedback in Equation (1) are 102 the size of the option value. These are unacceptable results that simply have no realistic interpretation. The nonlinearities of Equation (1) along with the non-continuous characteristics of the payoff function at
are introducing instabilities that can’t be handled by traditional FDM approaches. As explained above the traditional payoff function can produce instabilities in valuation estimates as can be seen in Figure 1 and Figure 2. [48] suggest using the payoff function
for puts (and
for call options).2
![]()
Figure 1. Impact of liquidity on put option value using finite differences.
![]()
Figure 2. Impact of liquidity on call option value using finite differences.
Figure 1 and Figure 2, despite the clear numerical limitations of FDM, still illustrate some interesting findings. The impact of the feedback parameter when
is greater than 1 is intuitive. As the parameter increases the pricing impact generally increases. If we interpret this parameter as a proxy for liquidity, then this result matches what we expect. That is, as illiquidity increases the impact on the options value should be more significant because the ability to create a hedge increases. As the term
approaches 1, the impact on pricing becomes more significant as well. The interplay between these two terms may be at the source of some of the numerical properties presented above.
Figure 1 illustrates the instability of the numerical approximation to a put option as a function of the liquidity parameter when using the finite difference approach to the Equation (2).
Figure 2 illustrates the instability of the numerical approximation to a call option as a function of the liquidity parameter when using the finite difference approach to the Equation (2).
7. PSM Results
Figure 3 and Figure 4 present the results using our PSM approach for Equation (1) under varying feedback constants less than 1 for each option type. PSM eliminates the erratic numerical results found using the FDM. This is an extremely important result. PSM is far more stable. While the general impacts increase with illiquidity, that impact is economically more defensible using PSM than FDM. Also, the impacts are smooth functions with peak impact near
where we have the maximum impact of the nonlinearities. For
we see a monotonic pricing impact of illiquidity. This can be seen in both Figure 3 for the put option as well as Figure 4 for the call option. Note also, that the boundary issues for the put option entirely disappear using PSM.
![]()
Figure 3. Impact of liquidity (<1) on put option value using PSM.
![]()
Figure 4. Impact of liquidity (<1) on call option value using PSM.
Figure 3 illustrates that numerical approximation to a put option for small liquidity parameters. Using PSM the approximation does not suffer any instabilities and illustrates that the liquidity parameter has the most impact on valuation around the strike price. Moreover, as the parameter gets larger, the impact on value increases.
Figure 4 illustrates that numerical approximation to a call option for small liquidity parameters. Using PSM the approximation does not suffer any instabilities and illustrates that the liquidity parameter has the most impact on valuation around the strike price. Moreover, as the parameter gets larger, the impact on value increases.
Figure 5 represents the results for the put options and Figure 6 for the call option for values of
. These results are somewhat counter intuitive but consistent with the findings of [54] except there they artificially avoid the counter-intuitive results by imposing arbitrary properties onto the solution. As
increases to very large values we see that the pricing impact changes sign. Note also, that large
also has larger price impact relative to
. This is intuitive and matches earlier studies. While the impact is positive, the effect of the parameter matches that of above; larger parameters have larger impact. Monotonicity in both areas is maintained, however, suggesting that the result is not a numerical flaw but a property of the pricing dynamic. Analyzing Equation (1), we observed that the product
is introducing larger impacts as it gets closer to 1. The feedback parameter and the gamma term (
) may be interacting to produce results that increase the price impact more than expected. As the gamma gets larger at the money, the need to hedge the risk-free portfolio becomes far more frequent, so that any illiquidity has a larger impact on valuation. Future study might focus on the details of this interaction using the more robust pricing approach introduced here.
Figure 5 illustrates that numerical approximation to a put option for large liquidity parameters. Using PSM the approximation does not suffer any instabilities and illustrates that the liquidity parameter has the most impact on valuation around the strike price. Interestingly, the impact of the liquidity parameter on valuation changes sign as it increases.
Figure 6 illustrates that numerical approximation to a call option for large liquidity parameters. Using PSM the approximation does not suffer any instabilities and illustrates that the liquidity parameter has the most impact on valuation around the strike price. Interestingly, the impact of the liquidity parameter on valuation changes sign as it increases.
3Presenting payoff graphs for
is not useful since for each
, the payoff profiles are indistinguishable from one another.
Figure 7 and Figure 8 present the price impacts for
using the payoff profile for each option type.3 These more traditional graphics are consistent with earlier studies.
![]()
Figure 5. Impact of liquidity (>1) on put option value using PSM.
![]()
Figure 6. Impact of liquidity (>1) on call option value using PSM.
![]()
Figure 7. Impact of liquidity (>1) on put option payoff function using PSM.
![]()
Figure 8. Impact of liquidity (>1) on call option payoff function using PSM.
Figure 7 illustrates the impact of the liquidity parameter over a wide spectrum of underlying asset prices. It shows that the PSM approach offers a very smooth approximation across the entire spectrum, as well as illustrating the impact of the parameter on the put option value. As the parameter increases so does its impact on value relative to the Black-Scholes solution.
Figure 8 illustrates the impact of the liquidity parameter over a wide spectrum of underlying asset prices. It shows that the PSM approach offers a very smooth approximation across the entire spectrum, as well as illustrating the impact of the parameter on the call option value. As the parameter increases so does its impact on value relative to the Black-Scholes solution.
8. Conclusion
We have introduced the PSM framework as a viable approach to numerically (or symbolically) analyze the NPDE’s frequently found in the contingent claim valuation literature. To illustrate the strengths of PSM we evaluated a commonly found NPDE where a feedback term is introduced to the dynamics of the riskless hedge portfolio. Many of the shortcomings in other numerical and analytic approaches are avoided when using the PSM approach. It more easily deals with the discreteness found around the strike price which frequently causes many other numerical frameworks to become unstable. PSM can produce a symbolic as well as numeric solution each of which can be analyzed in terms of the Greeks and sensitivity to the solution. We illustrated this impact and the corresponding improvements using PSM. For those interested in a more stable and robust numerical scheme, and one that offers a complete set of first and second order solutions across the entire numerical grid, PSM offers a viable and accurate alternative. Future research might extend the approach to more esoteric applications found throughout the literature and in practice.
Appendix
We present how to use Picard iteration in Maple to generate a symbolic polynomial approximation for the solution V to Equation (2). Two results from Maple are presented below.
To apply Picard iteration to Equation (2), one must first write it as the equivalent integral equation
Then one must rewrite this in polynomial form. We do this by using the auxiliary variables
and the chain rule to obtain a system of integral equations that can be integrated one at a time in order in Maple. The result from differentiating Equation (2) and algebra is
where,
The second-degree symbolic approximate solution for
produced by Maple is

Maple generated the fourth degree symbolic approximate for
with
but the output was too large to present but the coefficient of
is

To emphasize, note that Maple can be used to generate the Greeks symbolically and then these can be analyzed in a variety of ways. Of course, as pointed out in this paper and the linear paper, all the Greeks can be generated numerically using the discrete PSM J polynomials approximating
at
.