Non-Negativity Preserving Numerical Algorithms for Problems in Mathematical Finance

We give a study result to analyze a rather different, semi-analytical numerical algorithms based on splitting-step methods with their applications to mathematical finance. As certain subsistent numerical schemes may fail due to producing negative values for financial variables which require non-negativity preserving. These algorithms which we are analyzing preserve not only the non-negativity, but also the character of boundaries (natural, reflecting, absorbing, etc.). The derivatives of the CIR process and the Heston model are being extensively studied. Beyond plain vanilla European options, we creatively apply our splitting-step methods to a path-dependent option valuation. We compare our algorithms to a class of numerical schemes based on Euler discretization which are prevalent currently. The comparisons are given with respect to both accuracy and computational time for the European call option under the CIR model whereas with respect to convergence rate for the path-dependent option under the CIR model and the European call option under the Heston model.


Introduction
Stochastic differential equations (SDEs) are fundamental in mathematical finance. Particularly, they serve as models for describing the evolution of certain financial variables, such as the stock price, interest rates or volatility of an asset.
Although there are extensive studies on SDEs, explicit solutions of SDEs are How to cite this paper: Yuan, Y. (2018) Non-Negativity Preserving Numerical Algorithms for Problems in Mathematical Finance.

Problem Statement and Motivation
In some cases, the Itô-type SDEs of the form are well-defined only with certain boundary conditions. For example, the Cox-Ingersoll-Ross (CIR) [5] model and the Heston model [6] which were proposed to describe the short rate of interest and stock price dynamics respectively are well-defined in the domain X t = implies an absorbing or reflecting state. Precisely, the CIR process can be expressed in the where ( ) V W t is a Wiener process and , , κ θ σ are positive constants. The Heston model is a two-factor model of the form:  [7], if 2 d ≥ , the boundary ( ) 0 V t = is unattainable and the strong solution is strictly positive; if Second, as the square-root term in CIR process avoids the possibility of negative values, Euler-Maruyama [1] and higher order Taylor-type methods [1] [2] cannot be applied in this case as they lead to negative values of ( ) X t in practice. The numerical methods which preserve non-negativity are preferred. Actually various numerical schemes have been devised to meet this requirement, among which the schemes based on Euler discretization normally gives high efficiency, but their strong convergence rate and discretization errors are difficult to establish. In contrast to the Euler-type methods giving strong convergence of order 1 2 ε − , Moro and Schurz [9] provide the splitting-step methods of strong convergence at least of order 1. Although they are more costly with respect to computational time, they give higher accuracy and convergence rate. With this motivation, we would like to evaluate some important financial instruments using splitting-step methods with comparison to Euler-type numerical schemes.
We apply them to option pricing including a path-dependent option valuation comparing both accuracy and computational cost.
Third, exact simulation methods exist for both the CIR process (see Glasserman [10]) and the Heston model (see Broadie and Kaya [11] [12]). The drawback of these exact simulation methods is that they are very time-consuming.
They are competitive when one just simulates the process at one time (or few times), for example to price a European option with a Monte-Carlo algorithm.
On the contrary, they are too slow when one has to simulate the process along a time-grid, for example to compute a path-dependent European option. Thus, the splitting-step methods may have advantage on pricing path-dependent options as they are more efficient.

Review of Research Status
For the CIR process, a number of numerical schemes based on implicit time-stepping integrators have been devised for the case of unattainable boundary condition, see for example Alfonsi [13], Kahl and Schurz [14] and Dereich et al., 2012 [15]. There also are some splitting methods that serve for both boundary conditions without any restriction of parameters, such as Alfonsi [16], Moro and Schurz [9]. The latter paper gives a splitting-step method which high- The exact simulation method for CIR model exists, see Glasserman [10], but it requires more computational time than discretization schemes (up to a factor 10). This was analyzed in Alfonsi [13]. For the Heston

Summary of the Paper
In the remaining paper, we introduce the general structure and properties of splitting-step algorithm proposed by Moro and Shurtz in [9] in Section 2, including an efficient sample manner for generating a non-central chi-square distribution random number.
In Finally we conclude and point out the future work in Section 5.

Construction of the Algorithms
The splitting-step algorithm [9] has general structure as followed. Suppose that the more general equation to be integrated is of the form where ( ) W t is a standard Wiener process? Then decompose the above equation into two subsystems where it's required that we know the exact strong solution for ( ) 1 X t or the conditional probability We adapt the example in [9] to see how to obtain ( ) 1 X t through the conditional transition probability , then the conditional probability distribution is given by is the Dirac delta function and I is the modified Bessel function of the first kind with index ν which is given by  . Then we can sample the conditional probability distribution function using the rejection or inverse methods but it's computational expensive.
There is a more simple way to obtain ( ) and K is chosen from a Poisson distribution with mean 2 λ and i z are independent Gaussian random numbers with zero mean and unit variance.
Generally, a non-central chi-square distribution random variable with d degree of freedom and non-centrality parameter λ whose probability density function is given by see Johnson et al., 1994 [30]. This distribution is properly defined for d positive, and was extended to the case 0 d = by Siegel [31], afterwards was extended to the case for 2, 4, d = − −  by Moro and Schurz [9]. They used the fact that the distribution of (0.1) can be expressed as a mixture of central χ 2 variables with Poisson weights, then it can be sample as follow: Choose K from a Poisson distribution with mean = 2, then ( ) can be sampled using any standard random number generator of the 2 d χ distribution. This sampling is used specially when λ is small. For the case of large λ , other approximations (see e.g. Johnson et al., 1994 [30]) might be taken into account for improving the speed of the splitting-step algorithms.

Strong and Weak Convergence of Splitting-Step Methods
be any random partition of the given time interval [ ] 0,T with sufficiently small maximum step size Then the time discretized approximation X ∆  of a continuous-time process X, is said to be of general strong order of convergence γ to X at time T if there exists a positive constant C, which does not depend on ∆, and a 0 0 δ > such that the following strong error ( ) Along with the strong convergence, the weak convergence can be defined. A discrete-time approximation Y ∆ is said to converge with weak order 0 β > to X at time T as 0 ∆ → if for each smooth function g of polynomial growth there exists a constant g C , which does not depend on ∆ and The splitting-step algorithm has been shown that has strong and weak order 1.0 of convergence under certain assumption of the coefficient functions in [9].
We quote the convergence theorem here: Recall that the original SDE is We refer to the splitting Theorem Assume that the coefficient functions with exclusively uniformly bounded derivatives are such The proof and a general theorem on L 2 -convergence based on variation-of-constants formula (VOP) see Moro and Schurz [9]. We check the weak convergence with numerical experiment by taking g the identity function in section 3, where we present numerical studies of various strategies for integrating SDEs.

CIR Model and Its Properties
The famous Cox-Ingersoll-Ross (CIR) model was introduced in 1985 by John C.
Cox, Jonathan E. Ingersoll and Stephen A. Ross in [5] as an extension of the Vasicek model. It serves to describe the evolution of interest rates. This model specifies that the interest rate dynamics follows the stochastic differential equation, which is also called CIR Process: with k, θ and σ strictly positive parameters and ( ) W t a Wiener process. The parameter κ determines the speed of adjustment, θ is the long-run mean and σ is the so-called volatility to volatility. The drift factor, is exactly the same as in the Vasicek model. It ensures mean reversion of the interest rate towards the long-term value θ, with speed of adjustment governed by the strictly positive parameter k. According to Feller's boundary criteria, see Feller [32], different k, θ and σ values may produce distinct behavior at the boundary κθ σ ≥ , the upward drift is sufficiently large to make the boundary unattainable, i.e., the solution is always positive ( ) 0 The boundary becomes attainable, but it is strongly reflecting. That is, when a sample path reaches 0, then it returns immediately to the positive domain in a reflecting manner.
Since for all positive values of k and θ, the standard deviation factor ( ) V t σ rejects the possibility of negative interest rates, integration strategies must preserve non-negativity. Otherwise, it not only lacks any possible interpretation in the context of finance but also could induce severe errors in option valuation.

Simulation Schemes for the CIR model
We now turn to the simulation of CIR model (12). The exact simulation method for CIR model exists, see Glasserman [10], but it's very time-consuming. A simple Euler discretization of the CIR process may produce negative values of ( ) Practitioners often choose either a absorption or reflection fix whenever the process attains a negative value. That is to say, they deal with it by either setting [33]). Other direct approaches base on forced Euler-Maruyama approximations by slight modification of the model are also prevalent, see for example Higham and Mao [20], Deelstra and Delbaen [17], Lord et al., 2008 [21], which the full truncation method in Lord et al., 2008 [21] in practice has been proved to be the leading method of this class. Apart from Euler-type schemes, the splitting-step methods we introduced in the foregoing section can be applied to CIR model without any modification of the SDE. Moreover, the non-negativity is preserved during all the time interval.

Numerical Schemes Based on Euler Discretization
Here we outline the numerical schemes based on Euler discretization which we mentioned above. Let Although it can be applied to the simulation, it leads to negative values of ( ) V t in practice as we can see from Figure 1 • In Lord et al., 2008 [21], when they referred the method Higham and Mao, they had one more step which was taking the absolute value of the right side of (3.2), we call this method Higham and Mao complemented, which is as follow • Deelstra and Delbaen [17] proposed an Euler-type scheme by taking positive part of the component inside the squared-root, this scheme is called partial truncation in Lord et al., 2008 [21] and has the form as follow • Full truncation was devised by Lord, Koekkoek & Van Dijk in [21] which is the leading method of these classes. The difference between full truncation and partial truncation is the treatment of the drift term, where the full truncation method has one more x and can be expressed as follow We refer the above four numerical schemes as Higham and Mao, Higham and Mao complemented, partial truncation and full truncation respectively.

Splitting-Step Algorithms Applied to CIR Process
According to the general structure of splitting-step algorithm in Section 2, we split the SDE (3.1) into two equations as follows.
We note that the process defined by (3.13) is a κθ-dimensional squared Bessel process (BESQ) with index of the process = in the case of a regular boundary, where I ν is the modified Bessel function of the first kind with index ν. Then we notice that the transition density  can be written in terms of a non-central 2 χ distribution.
Recalling the probability density function (PDF) of a non-central chi-square distribution random variable with d degree of freedom and non-centrality parameter λ:   V t t + ∆ has the initial data of (3.14) and integrate it by any deterministic numerical method of at least order 1. In our project, we practically sample the ( ) In Figure 1, three trajectories for the CIR process are shown with a comparison to Higham-Mao which is simply taking absolute value of ( ) V t in the squart-root term. Here we choose 1, 1, 2 κ θ σ = = = , so that the boundary is attainable. The figure depicts that the splitting-step algorithm preserves non-negativity for positive initial data while Higham-Mao produces negative values occasionally. Thus splitting method is preferable as the CIR process is not defined for negative values.

Weak Convergence
In this section we verify that the splitting-step algorithm has weak convergence of order 1.0 on the interval [0,T] for a fixed finite, nonrandom terminal time 0 T > .
We take g the identity function and use the model Then split it in this way: Taking the expected value of both sides of (3.19) reads: Notice that the second term on the right side Then we apply the splitting-step method to obtain the numerical approxima- Figure 2 depicts the error Versus decreasing uniform step size ∆t, which shows that the method has weak order 1.0 while using constant step sizes ∆t. Moreover, we compare to Higham-Mao complemented, partial truncation and full truncation for the same equation. Figure 2 shows that splitting methods has the highest precision among these four algorithms.

Option Valuations
Among the variety of financial derivatives, the option is one of the most important financial instruments. In current financial markets, there are mainly four kinds of options: American option, European option, Asian option, and Barrier option. In our work, we only focus on pricing European options, which can only be exercised on the maturity date whereas an American option can be exercised at any time before expiration.

A European Call Option
A European call option on an asset is a contract that allows the buyer to buy (but not obligate) this asset at the price K at time T. The number 0 K > is called the exercise (or strike) price and 0 T > the maturity (or expiration) time. We pay K to buy the asset at time T if ( ) as we can buy it cheaper than K. Thus at time T, the call option has the value Suppose that C T is the numerical approximation of C T , then the error reveals the precision of an algorithm. Since the exact solution C T no explicitly known, we compute a numerical reference solution with a very small t and regard it as C T . Here we use the numerical approximation of C T with 4 10 t = computed by splitting method as reference. Figure 3(a) depicts the comparison of the four algorithms with respect to the precision of calculating call option for each step size t and Figure 3(b) shows the computational cost versus precision among them. We observe that Higham-Mao complemented, partial truncation and full truncation work better on this call option pricing as they have more precision and higher speed.

A Path-Dependent Option
Apart from the typical European call option in the previous section, more exotic options exist in the market (see, for example [35] [36]). Here we consider one of path-dependent options: up-and-out call option (hereafter referred as U&O option). An up-and-out call option pays off the usual at expiry unless at any time during the life of the option the underlying asset reaches certain level (from below, obviously) then it is said to knock out and become worthless. The Figure 4 explains the payoff of each underlying asset. We assume that a underlying will be knocked-out once it reaches L from below at any time prior to expiry time T. Then the underlying asset corresponds to the blue one has payoff 0 as it's knocked-out.
As we say that an U&O option expires worthless if the asset price touches some barrier L from below, say, at any time prior to expiry, where L is larger than the present asset value 0 V , the price of U&O call option at expiry time T is given by  Let t u ∆  which depends on the time step ∆t be a numerical approximation of an exact value u. The numerical method said to be of order p means that there exists a number C independent of ∆t such that Hence, we can get an estimate of the order of accuracy p after computing The relative differences indicate the convergence rate of an algorithm. The faster converges a scheme (has higher p value), the faster relative differences reduce to 0 as ∆t tends to 0. Now we define the relative difference of U&O call option value as:  computed by splitting by the fit of existing data. This has implied that splitting convergent the fastest as it's the first one to reach the SEM level. Figure 5 . But along the trend of smaller relative differences, splitting performs better and better and costs the least since certain value of relative difference. For example, with the same parameter setting, we can deduce that splitting will cost the least for the cases of & 0.001 from the tendency of these four algorithms.

Heston Model Basics
The Heston model [6] is defined by the coupled two-dimensional stochastic dif- The parameters μ is the rate of return of the asset. θ is the long variance, which means as t tends to infinity, the expected value of ( ) V t tends to θ. κ is the rate at which ( ) V t reverts to θ σ is the so-called vol of vol, which determines the variance of ( ) V t . We note that the variance (4.2) follows a CIR process.
For ( ) V t it has the same boundary behavior as we mentioned in §3.1.
There is an additional condition for ( ) S t which is iii) ( ) S t has an absorbing barrier at 0.

Path Simulation
There are plenty of methods can be used to simulate Heston model. Broadie and Kaya [11] [12] derived an exact simulation method without bias, but highly time-consuming. This is analyzed in Broadie and Kaya [21], and Lord et al., 2008 [21]. The exact scheme is competitive when one simulates the process just at one time (or few times), for example to price European options with a Monte-Carlo algorithm. On the contrary, they are drastically too slow for simulating the process along a time-grid, which occurs when computing path-dependent options prices. Kahl and Jäckel [22]

A European Call Option with Heston Model
Closed form semi-analytical formulae for plain vanilla option prices have been derived in [6]. However, these formulae require the evaluation of logarithms with complex arguments during the involved inverse Fourier integration step.  Obviously the splitting-step method has the highest order of accuracy p. Although the splitting method is more time-consuming than the others for each ∆t, and with respect to relative differences, it costs the most at the beginning as for the case 2 Hes t ∆ >  illustrated in Figure 6(b). But it costs the least computational time if one restricts that the relative differences should be less than a certain value. For example, if one requires that the relative difference should be less than 0.10 with the same parameter setting, splitting practically works much faster than the others. The Figure 6(b) gives the evidence.
More precisely, we analyze the computational cost by taking the case 0.065

Conclusions
A number of existing numerical algorithms for integrating SDEs cannot be used to simulate certain financial models which require non-negativity, such as the CIR model, the Heston model that we have seen. A slight modification of Euler-Maruyama [1] method conducts to various Euler-type numerical schemes.
They are simple to implement and efficient enough, but have lower accuracy compared to exact simulation schemes. Sometimes they also give (4.5). The vertical dash line corresponds to 0.065 Less precision is compared to semi-analytical numerical schemes such as the splitting-step methods which we are analyzing in our project. Since the exact simulation schemes are thoroughly too slow for simulating a process along a time-grid, splitting-step methods preponderate them with respect to efficiency.  The splitting algorithms heavily rely on the exploitation of the specific structure of original system (2.1). One can decompose the original system into two subsystems appropriately for which either one knows the explicit solution or the conditional transition probability of one of the subsystems and integrate the remaining one by deterministic numerical methods of at least order 1. In this way, it preserves the non-negativity and a maximum of convergence order 1.0 both in strong and weak sense. For the analytical solution of SDEs, an extensive list of known one can be found in textbooks such as Kloeden and Platen [1].

Future Work
Since derivatives are one of the three main categories of financial instruments where the other two being stocks (i.e., equities or shares) and debt (i.e., bonds and mortgages), and as we said that the splitting methods are much more efficient but has less accuracy than exact simulation methods which are not suitable for pathwise simulations, then one direction of our interest is applying splitting-step method to other path-dependent options, probably with comparison to other efficient algorithms like what we have done in §3.3.2.
Another direction is the application of splitting-step models to portfolio problems which are aiming to find the optimal investment strategy of an investor. Following the optimal portfolio strategy leads to the maximum expected utility of the terminal wealth. For example the classical Merton's problem [37] and mixed stock & bond portfolio problem like in [38] are under consideration.
As Moro and Schurz [9] have mentioned that the splitting-step scheme is applicable to multi-dimensional problems, we intend to encounter some suitable splitting manners to deal with them which in particular can be SDEs with linear, decoupled, commutative or diagonal noise terms. Furthermore, non-linear partial differential equations, for instance, Super-Brownian motion interest us to verify its properties by implement splitting-step methods.