A Cubic B-Spline Collocation Method for Barrier Options under the CEV Model

In this paper, we construct a new numerical algorithm for the partial differential equation of up-and-out put barrier options under the CEV model. In this method, we use the Crank-Nicolson scheme to discrete temporal variables and the cubic B-spline collocation method to discrete spatial variables. The method is stable and has second-order convergence for both time and space variables. The convergence analysis of the proposed method is discussed in detail. Finally, numerical examples verify the stability and accuracy of the method.


Introduction
With the further acceleration of economic globalization, the price turbulence in the financial market is becoming more and more serious, and various risks caused by it have become a hot issue concerned and studied by economists and investors.In recent years, the derivatives market has become increasingly important in the financial and investment fields.Options have become an important tool for companies and investors to hedge assets and investment portfolios to control the risks brought by stock price fluctuations.There are many types of options, such as European options, American options, Asian options, barrier options, exotic options, lookback options, etc.European options are options that one party of the call option must exercise on the expiration date of the option.American options and barrier options are path-dependent options, and the return of these two options depends on the volatility path of the underlying asset price, not just the asset price on the maturity date.
The famous Black-Scholes option pricing model was proposed by Fischer Black and Myron Scholes in 1973 [1].This model assumes that the underlying asset price follows the geometric Brownian motion, i.e., its distribution is logarithmic normal, and when it focuses on option pricing under the premise that its volatility is constant, it can meet option pricing in most cases.Since the required preconditions are too strict, the conditions proposed by the Black-Scholes model do not conform to the reality.In 1975, Cox [2] first considered the characteristics of the relationship between the volatility of assets winning the bid in the market and their price level, and established the Constant Elasticity of Variance (CEV) model.Schroder's [3] empirical research showed that the CEV model can more accurately describe the changes in the price of underlying assets.Therefore, applying the CEV model to option pricing can obtain a better solution.
The option pricing model used in this section is the CEV model.The stock price process satisfies: where δ is a constant coefficient, 0 δ 1, µ is expected rate of return, W t is standard Brownian motion.
The differential equation that the option price satisfies is where U (S, t) is the value of options at the asset price S and at time t, T is the maturity date, σ is the volatility and r is the risk-free interest rate.Barrier options are a kind of European option, and their final benefit depends not only on the price of the underlying asset on the expiration date of the option, but also on whether the price of the underlying asset reaches a certain level (barrier B) during the entire option validity period.When the price of the underlying asset reaches the barrier, the option terminates.This type of option is called a knock-out option.If the underlying price of the original asset is higher than the barrier, it is called a down knock-out barrier option; if the underlying price of the original asset is lower than the barrier, it is called an up knock-out barrier option.When the price of the underlying asset reaches the barrier period, the option becomes effective, and this type of option is called a knock-in option.Similarly, there are down knock-in options and up knock-in options.European options can be divided into call options and put options.
Let us consider the up-and-out put barrier options, the initial condition and boundary conditions where K is the exercise price.Due to the fact that the initial condition is not smooth, the obtained solution is not smooth enough to perform the convergence of numerical approximations.We replace max(K − S, 0) by a smooth function ψ(S).We obtain the initial condition and boundary conditions Using the transformation: S = exp(x).We obtain with initial condition and boundary conditions where Now, to solve the problem given in ( 9)-( 12), we truncate the infinite domain (−∞, x max ) into a finite domain (x min , x max ), where x min is a sufficiently small negative real number.Then, we consider a general form of ( 9)-( 12) on the truncated region Ω = (x min , x max ) × (0, T): where with initial and boundary conditions In order to ensure the existence of a unique solution v ∈ C( Ω) ∩ C 2,1 (Ω) of the above problem, we assume that [4] ∂ l+n v(x, t) ∂x l ∂t n K on Ω, 0 n 3 and 0 l 4.
There are numerous studies about the numerical methods of solving barrier options under the CEV model.For example, Thakoor et al. [5] provided a highly accurate and efficient numerical scheme for barrier options under CEV.In order to achieve the fourth order convergence, it is necessary to refine the grid at both the strike and barriers, and obtain oscillation-free hedging parameters; a full implicit time integration technique with extrapolation is employed.Ahmadian and Ballestra [6] developed a numerical method to price discrete barrier options based on an underlying described by the constant elasticity of variance model with jump-diffusion (CEVJD).This method is allowed to approximate differential terms and jump integral by means of two different numerical techniques.Dias [7] has developed two new methods for pricing and hedging European barrier option contracts under the Jump to Default Extended Constant Elasticity of Variance (JDCEV) model, i.e., the stopping time method based on the first passage time densities of the underlying asset price process through the barrier level, as well as a static hedging portfolio method in which the barrier option is replicated by a portfolio of plain-vanilla and binary options.Guo and Chang [8] proposed an accelerated static replication method for continuous European barrier options by using the repeated Richardson extrapolation technique for Romberg sequences.
The main motivation of this work is to develop a third-order numerical method for solving the CEV model given in ( 13)- (16).In this method, the Crank-Nicolson scheme is used to discrete the temporal variables, and the cubic B-spline configuration method is used to discrete the spatial variables.The stability of this method is analyzed.In addition, we prove that the proposed method is second order convergent in space and time.The the-oretical results are verified by numerical experiments and the applicability of the method is illustrated.B-spline was originally proposed by Schoenherg [9] in 1946.B-spline has strong local control capabilities, as well as many advantages such as convexity preservation, convex hull, geometric invariance, and variability reduction.Nowadays, B-spline has become one of the most important numerical methods used to study and solve partial differential equations.For example, Shallu et al. [10] proposed an improved extrapolation collocation method based on cubic B-splines.The method is applied to a class of self adjoint singularly perturbed boundary value problems, with the equispaced discretization of the domain.Since the standard B-spline collocation method cannot achieve the optimal convergence order, the cubic B-spline interpolant and its higher-order derivatives are corrected later.The proposed technology can be applied to high-order linear, nonlinear problems, and partial differential equations.Shallu and V.K. Kukreja [11] obtained a numerical solution to the nonlinear modified Burgers equation using an improvised collocation technique with cubic B-spline as basis functions.In this method, the cubic B-spline is forced to satisfy interpolant conditions with some specific end conditions.Shallu and V.K. Kukreja [12] uses an improvised cubic B-spline collocation technique to solve the Benjamin-Bona-Mahony-Burgers (BBMB) equation due to the inability to achieve optimal accuracy and convergence order using the standard B-spline collocation method.Therefore, improvised cubic B-splines are formed by making a posteriori corrections to cubic B-spline interpolants and their higher-order derivatives.We note that numerical methods based on B-spline basis functions have been used to solve option pricing problems.For example, the authors of [13][14][15] used the B-spline basis function for solving the generalized Black Scholes equation.
The rest of this paper is organized as follows: In Section 2, the time variable is discretized by the Crank-Nicolson format, and the convergence analysis is carried out.In the Section 3, we describe spatial discretization by the cubic B-spline collocation method, and discuss the stability and convergence analysis of this method.Numerical experiments are carried out in Section 4. The results will be discussed in Section 5.

Temporal Discretization
In order to solve the CEV model given in ( 13)-( 16) more conveniently, we consider a uniform partition of the domain Ω = (x min , x max ) × (0, T) with node points (x l , t n ), where We first discretize the problem ( 13)-( 16) in the temporal direction.

Description of the Crank-Nicolson Method of Time Variable
We use the Crank-Nicolson method to discretize Equations ( 13)-( 16).This discretization yields the following system of differential equations: Here, v n (x) approximate the exact solution v(x, t n ) at the nth time level t n and

Convergence Analysis
Next, we will conduct convergence analysis for the case of using the CN model discrete temporal variables.
The local truncation error of the temporal discretization method ( 17) is defined by where ṽn+1 is the approximate solution of the problem ( 17)-(20) obtained after one step of the semi-discrete scheme taking the exact value v(x, t n ), instead of v n as the initial data, i.e., We define the global truncation error E n+1 at time t n+1 as x (w) is the solution 'u' of: , where R( L x ) ⊆ C(Ω x ) denotes the range of the linear operator L x .

Lemma 2.
If then there holds e n ∞ Kτ 3 (29) where K is a constant independent of N.
Proof.Because |v ttt | K and v is a smooth function, hence by Taylor series expansion, we obtain Using the above extension, we obtain Using relations ( 13) and ( 21) gives By relations (31) and (32), we obtain Now, because This means It follows by (33) that By using relation (24), we obtain By directly application of Lemma 1, we obtain In addition, we can decompose the global truncation error at time t n in the following form By relations ( 17) and ( 24), we obtain (34) and (35) gives the recurrence relation as Here, we use the notation Now, according to Gonzalez and Palencia (see, [17]) for the operator 'R n ≡ R(τL n−1 x )', where R(z) is a rational A-acceptable function (see, [18]), such that Thus, we obtain This imples that (30) follows.
We can see that the discrete scheme is consistent by Lemma 2. Thus we obtain the following result with the help of Lax equivalence theorem (see, [19]) and Lemma 2.

Description of Cubic B-Spline Method
First, we consider a uniform partition π : over the finite interval Ωx = (x min , x max ).Then we extend the uniform partition π in order to construct the cubic B-spline for the partition Let Q n+1 (x) be the approximation to the exact solution v n+1 (x) at the (n + 1)-th time level.Then, Q n+1 (x) can be defined by where c n+1 k are unknown parameters.Following [20], we can obtain the representations of cubic B-splines B k with the required attributes: Each basis function is twice continuously differentiable on the entire real line.
If we consider B 3 (π By using (37), the values of basis functions B k (x), B k (x), B k (x) at the grid points x l are listed in Table 1.By using (36) and Table 1 xx at the mesh points x l are obtained as Table 1.The values of basis functions B k (x), B k (x), B k (x) at the grid points x l .

Formulation of Cubic B-Spline Collocation Method
In this part, we construct the cubic B-spline collocation method for the solution of the problem ( 17)-(20).Firstly, we write Equation ( 17) as: where the boundary conditions Then, we force the function Q n+1 (x) to satisfy Equations ( 41) and (43) at the node points x l , according to the cubic B-spline collocation technique.Thus, we have Substituting the expressions for Q n+1 (x l ), Q n+1 x (x l ) and Q n+1 xx (x l ) from Equations ( 38)-(40), respectively into Equation (44), we obtain a system of m + 1 equations in m + 1 unknowns where By using (45), we obtain By using Equations ( 47)-(48), we eliminate c n+1 −1 , c n −1 , c n+1 M+1 and c n M+1 , such that where The system produced by Equations ( 49)-( 51) can be expressed in matrix form as follows where One can easily see that for sufficiently large m, the matrix P is strictly diagonally dominant by rows.Hence P is invertible (see, [21]), and We need obtain the initial vector T to solve the system of matrix (52).By solving the following interpolation problem can find the vector C 0 where Thus we obtain the following system of matrix where , . Now, since matrix A is diagonally dominant with strict diagonal dominance in all rows except the first row and the last row, it can be concluded from Gershgorin's theorem (see [22]) that A is nonsingular.Therefore, the vector C 0 giving the initial vector C 0 can be uniquely found.

Stability Analysis
In this subsection, we use the Von-Neumann stability method to analyze the stability of the scheme defined by (52).
Proof.We set f n+ 1 2 (x) = 0 to establish the stability result of method (52).We define c n l as where i = √ −1, ρ is the time related parameter and the exponential term represents spatial correlation.Substituting Equation (53) into Equation (46), we obtain µ n+1 e ilρh e −iρh p + q + e iρh r = µ n e ilρh e −iρh − 2 Re-write above equation as For the convenience of calculation, we introduce E 1 , E 2 and E 3 as By reducing Equation (55), we obtain with By taking the absolute value on both sides of Equation (57) we obtain Thus, we can see that the scheme is unconditionally stable by Equation (59).

Convergence Analysis
Lemma 3. Let S 3 (π) := set of all functions s(x) ∈ C 2 (Ω x ) that reduces to cubic polynomials on each subintervals (x l , x l+1 ), and, there exists a unique cubic spline s(x) ∈ S 3 (π) solving the interpolation problem- where, v n (x) is the solution of the discretized problem (41)-( 43) at the time level n.The function 's' is called the cubic spline interpolate to v n (x).
Theorem 3. The collocation approximation Q n (x) from the space B 3 (π) to the solution v n (x) of the boundary value problem (41)-(43) satisfies the following error estimate for sufficiently small h (i.e., for sufficiently large m) where K is a positive constant.
Proof.In order to find error estimate, let Q n (x) be the unique spline interpolate from Q 3 (π) to the solution v n (x) of the boundary value problem (41)-(43) at the nth-time level. Here, The boundary value problem ( 41)-( 43) can be rewritten as where Therefore, according to Hall and Meyer error estimate (see, [24]), we can obtain where, Thus, ∀l, 0 l M, and, ∀n, 0 n N (68) At the (n + 1)th-time level, we write where g n+1 is an error function which order of magnitude is O(h 2 ) and ψ n+1 L , ψ n+1 R are error functions which order of magnitude are O(h 4 ).
Since v n+1 (x l ) = Q n+1 (x l ); ∀l, 0 l m, therefore we can re-write the system (70)-(72) in the following matrix form where, By using ( 52) and (73), we obtain A simple use of Equation (74) and bound of P −1 ∞ derived in Section 3.2 gives where In addition, with the help of Equations ( 71), ( 72) and (75), we obtain Now, By using B k (x l ) given in Section 3.1, we can easily see that In the following triangular inequality, we discuss the relations (66) and ( 78) Theorem 4. Let v(x, t) be the solution of the initial boundary value problem ( 13)-( 16) and Q n (x) be the collocation approximation to the solution v n (x) after the temporal discretization stage.Then under the hypothesis of Lemma 2 and Theorem 3 the error estimate of the totally discrete scheme is given by v where K are positive constants.

Numerical Experiments
In this section, we mainly use RStudio version 4.2.0 to verify the feasibility of the obtained pricing model (52).
Consider the following issue of PDE of up-and-out put barrier options under the CEV model: subject to the initial and boundary conditions Figure 1 shows the change trend of Spatial Variable x and Option Price Q for δ = 1/4, 1/3, 1/2 when r = 0.05, σ = 0.05.When δ takes different values, it is noted that the changes in option prices are basically consistent, indicating that the obtained pricing model is feasible.Figure 2 shows the change trend of Spatial Variable x and Option Price Q for r = 0.05, 0.075, 0.1 when δ = 1/2, σ = 0.05.When the value of x approaches 2, the option price gradually increases as r increases, and the change trend of Spatial Variable x and Option Price Q is roughly the same in the rest of the time.Figure 3 shows the change trend of Spatial Variable x and Option Price Q for σ = 0.05, 0.075, 0.1 when δ = 1/2, r = 0.05.When the value of x is less than −1, the option price gradually increases as σ increases, and the change trend of Spatial Variable x and Option Price Q is roughly the same in the rest of the time.Overall, the option price gradually increases with the increase of spatial variable x, and the rate of growth also gradually increases.Figure 4 shows the change trend of Time Variable t and Option Price Q for r = 0.05, 0.075, 0.1 when δ = 1/2, σ = 0.05.The option price gradually increases as the time variable t increases.As can be seen from the figure, as r increases, the slope of the line gradually steepens, and the option price increases more quickly.When δ = 1/2, σ = 0.05 and r = 0.05, the surface plot of the approximate solution for M = N = 30 is given in Figure 5.   Based on the studied option pricing model, in this section we considers the relationship between the time variable x and the spatial variable t and the option price, and select different parameters δ, r and σ.Based on the selected parameters, the trend chart of the option price was obtained, and the changes in the image were observed.The numerical simulation results were consistent with theoretical analysis.

Conclusions and Dicussion
In this paper, we study the cubic B-spline method for the barrier option problem under the CEV model.The time variable is discretized using the Crank-Nicolson method.On this basis, the spatial variable is discretized using the cubic B-spline collocation method to obtain a pricing model for barrier options.The time discretization of barrier options is analyzed for convergence, and the space discretization is analyzed for stability and convergence.The trend of log price of stock and option price for each parameter is analyzed in RStudio version 4.2.0, and the results show that cubic B-spline is feasible for solving barrier option problems under CEV models.
At present, there are still some shortcomings in this article.When truncating the infinite domain into the finite domain, some errors may occur.However, according to references [13][14][15], these errors are controllable.Thus, in the data experiment section, we make the error f (x, t) = 0. Additionally, in this article, there is no comparison between the B-spline method used in this article and other methods for option pricing, which is also a limitation of this article.