On the Convergence of a Crank–Nicolson Fitted Finite Volume Method for Pricing American Bond Options

)is paper develops and analyses a Crank–Nicolson fitted finite volumemethod to price American options on a zero-coupon bond under the Cox–Ingersoll–Ross (CIR) model governed by a partial differential complementarity problem (PDCP). Based on a penalty approach, the PDCP results in a nonlinear partial differential equation (PDE).We then apply a fitted finite volumemethod for the spatial discretization along with a Crank–Nicolson time-stepping scheme for the PDE, which results in a nonlinear algebraic equation. We show that this scheme is consistent, stable, and monotone, and hence, the convergence of the numerical solution to the viscosity solution of the continuous problem is guaranteed. To solve the system of nonlinear equations effectively, an iterative algorithm is established and its convergence is proved. Numerical experiments are presented to demonstrate the accuracy, efficiency, and robustness of the new numerical method.


Introduction
Pricing and hedging of interest rate derivatives, such as bond options, have been increasingly attracting much attention by mathematicians and financial engineers in the last two decades [1][2][3]. Usually, bond option has a bond as its underlying asset whose price depends on both interest rate and time. Compared to stock derivatives, the valuation of these types of options is more complicated and poses greater challenges. In this paper, we focus on pricing American options on a zero-coupon bond under the CIR model governed by a PDCP (1). Due to the early exercise feature, this PDCP is, in general, not analytically solvable, and then, the efficient numerical approximation methods are normally required [4].
In the past few years, the numerical solution of the American bond option pricing problem has been investigated by various approximation techniques, e.g., the lattice method, the finite difference method, the finite element method, and the fitted finite volume method. For more details, we refer the readers to [5][6][7][8][9][10][11][12][13][14] and references therein. However, it should be noted that due to the degeneracy of the nonlinear PDE, many research studies, e.g., [5,11,15], have been pointed out that most of the above methods are difficult to solve the PDE accurately and/or only convergent for certain combinations of parameters. Recently, the fitted finite volume method for pricing stock options governed by standard Black-Scholes equations has attracted much attention. is method was first used to price European stock options by Wang [16] and then generalized to other types of options by many researchers (see, e.g., [14,[17][18][19][20][21] and references therein). Actually, the fitted finite volume method is based on a finite volume formulation (cf. [22,23]) coupled with a popular exponentially fitting approximation technique. It has been shown that this method has greater success in pricing stock options and is gaining popularity because it can overcome the difficulty caused by driftdominated phenomena very well. Besides, it is well known that based on the penalty approach, the nonlinear penalized PDEs resulting from American options models are usually degenerated and drift-dominated, and hence, the fitted finite volume method is a natural way to overcome these difficulties [17]. e fitted finite volume method in [11,14] evaluates the American and European options on a zero-coupon bond under the CIR model, respectively, and the corresponding PDEs are treated by a fully implicit (i.e., backward Euler) time-stepping scheme. We note that the scheme has the firstorder convergence rate. To the authors' knowledge, based on the fitted finite volume method discretization, there is no theoretical result published about the Crank-Nicolson timestepping scheme which will provide the second order of convergence to price this problem for the time being. On this basis, in this paper, we derive a novel fitted finite volume method combined with the Crank-Nicolson scheme for pricing American options on a zero-coupon bond under the CIR model. Based on a penalty approach, the PDCP (1) results in a nonlinear penalized PDE, and then, we apply the fitted finite volume method in space, coupled with the Crank-Nicolson scheme in time. We show that this numerical scheme is consistent, stable, and monotone, and hence, the convergence of the numerical solution to the viscosity solution of the continuous problem is guaranteed. Further, to solve the nonlinear system effectively, we design an iterative method and prove that the method is convergent. Finally, some numerical experiments are implemented to verify the accuracy, efficiency, and robustness of our numerical method. e numerical results further show that the Crank-Nicolson fitted finite volume scheme provides the rate of convergence is of approximately second order. e paper is organized as follows: in the next section, we introduce American options on a zero-coupon bond model and give a penalty approach to the PDCP. en, a Crank-Nicolson fitted finite volume scheme is constructed in Section 3. In Section 4, the convergence of the numerical scheme is investigated. In Section 5, we propose an iterative algorithm for the nonlinear system and establish a convergence theory. Numerical tests are performed in Section 6. e paper ends with some conclusions.

Mathematical Model and the Penalty Approach
In this work, we assume the short-term interest rate (denoted by "r") structure is governed by the CIR model, i.e., r is governed by the following mean-reverting version of the square-root process: where dW is the increment of a Wiener process, θ is the long-term level of the short rate, κ > 0 stands for the reversion speed, and σ 2 r is the variance with σ > 0. In [24], it has been shown that the price P(r, t, s) of a pure discount bond with face value $1 at its maturity date s is given as follows: where and ζ is the market risk premium. At the maturity date s, the price of a pure discount bond is its face value, i.e., Now, let V(r, t) be the value of an American put option on a zero-coupon bond with striking price K, where the holder can receive a given payoff Λ(r, t) at the expiry date T. Introducing a time-reverse transformation τ � T − t, then the value V(r, τ) can be governed by the following PDCP (cf. [12]): a.e. in (0, +∞) × (0, T), where and the initial condition is given by and the boundary conditions are as follows: For computational purpose, it is necessary to restrict the short-term interest rate r in a finite region I � [0, r max ], where r max is a sufficiently large number to ensure the accuracy. us, the boundary conditions (8) become By using the penalty method (cf. [11,20,25]), we can transform the PDCP (1) into the following nonlinear penalized PDE: with the given initial and boundary conditions: where c is the penalty parameter satisfying c > 1 and [x] + � max 0, x { } for any x. It is well known that, by adding the penalty term, the above penalty approach is to force the positive parts [Λ − V c ] + in (10) to be close to zero as c becomes sufficiently large. Hence, the complementarity condition in (5) is approximately satisfied.

Remark 1.
Since the degeneracy of the diffusion operators L, the nonsmoothness of payoff functions, and nonlinearity of max ·, · { }, both problems (5) and (10) have no classical solutions, in general. erefore, we seek viscosity solutions to the PDCP (1) and penalized PDE (5). Detailed discussions concerning the existence and uniqueness of viscosity solutions to PDE in the context of financial mathematics are given in [26,27], which are defined as follows.
Definition 1 (continuous viscosity solutions [27]). u ∈ C(Ω) is a viscosity solution of where H is a continuous function satisfying the ellipticity condition: for any x ∈ Ω, Ω is a domain in R N , u ∈ R, p ∈ R N , if and only if ∀φ ∈ C 2 (Ω), if x 0 ∈ Ω is a local maximum point of u − φ, one has and ∀φ ∈ C 2 (Ω), if x 0 ∈ Ω is a local minimum point of u − φ, one has We note that, according to [11], the solution V c to penalized PDE (5) is an approximation to solution V of PDCP (1) in the viscosity sense. e convergence of V c to V as c ⟶ ∞ has been discussed by many researchers, e.g., Li and Wang [28] and Wang et al. [29].
erefore, we will concentrate on the numerical approximation of the viscosity solution to PDE (5). Besides, to avoid an overloading of symbols, we omit the subscript c of V c used in the above, but bear in mind that this V is the solution of penalized PDE (5).

The Crank-Nicolson Fitted Finite Volume Method
In this section, we will present a Crank-Nicolson fitted finite volume method for the nonlinear penalized PDE (5).

Semidiscretization Scheme.
Before proceeding to the Crank-Nicolson fitted finite volume discretization scheme, we first transform PDE (5) into the following conservative form: where e fitted finite volume method is based on the selfadjoint form (16). Now we first define two spatial partitions of I � [0, r max ]. Let I be divided into M subintervals: Applying the one-point quadrature rule to all the terms in (19) except the first term in the right-hand side, we obtain where l i � r i+1/2 − r i− 1/2 is the length of interval I * i and ρ(V) is the weighted flux density associated with V defined as follows: where V i and Λ i denote the approximation of V(r i , τ) and Λ(r i , τ), respectively. As the discretization of (20) and (21) is identical to that in [11,14], so we omit the details and only give the final discrete form here, i.e., the semidiscretization of the PDE (6) is given by and S is a (M − 1) × (M − 1) semidiscretization matrix given by where and V 0 and V M being equal to the given boundary conditions in (9), and η i � b i+1/2 /a.

e Crank-Nicolson Scheme.
We now consider the timestepping scheme to (22). Let τ n � nΔτ(n � 0, 1, . . . , N) be a set of partition points in [0, T] with subinterval length Δτ � T/N, where N is a positive integer. It has been shown in [11,14] that the fully implicit fitted finite volume scheme is only first-order accurate for pricing both European and American zero-coupon bond options under the CIR model. To improve the accuracy in time step, the second-order Crank-Nicholson scheme should be employed here. us, the Crank-Nicolson scheme of (22) can be written in the following matrix form: for where I is a (M − 1) × (M − 1) unit matrix and V n i denotes the approximations of V(r i , τ n ) at τ � τ n for i � 1, 2, . . . , M − 1.
We comment that in (26), the Dirichlet boundary conditions (9) have been incorporated, and the initial condition is also incorporated as the payoff function. Furthermore, from [11], we can know that the matrix − S is an M-matrix. So, combing the definition of P(V n+1 ), we further have the following result for the numerical scheme (26). [30], p. 91).

Remark 3.
eorem 1 implies that the Crank-Nicolson scheme (26) satisfies the discrete maximum principle, which guarantees that the discrete arbitrage inequality holds in option pricing theory.

Convergence of the Numerical Scheme
In this section, we investigate the convergence property of numerical scheme (26). First, the following lemma shows that scheme (26) is consistent.
Proof. From the fitted finite volume discretization, we can see that the consistency of the numerical scheme (26) relies on the consistency of ρ(V). Let u be a sufficiently smooth function and u h be the discrete approximation of u. It has been shown in [11,16] that the exact and the discrete flux yield where ρ h (V) satisfies for i � 0, 1, . . . , M − 1, and h � max max i h i , Δτ be the mesh parameter. Hence, it is easy to see that the consistency of the numerical scheme (26) is the consequent result. □ e stability result for the numerical scheme (11) is given as follows.

Lemma 2.
e numerical scheme (26) is stable as h ⟶ 0, i.e., Proof. Writing out (26) in component form gives It follows from β i ≥ 0, w i ≥ 0, χ i ≤ 0 that, when h is sufficiently small en, we obtain the following estimate: where − χ l − β l − w l � r l > 0. Combining (7) and (34) gives for all feasible n. Hence, the numerical scheme (26) is stable.
□ e monotonicity of the numerical scheme (26) is given in the following lemma.
i . Hence, the numerical scheme (26) is unconditionally monotone. □ It has been shown in [27] that the solution of the numerical scheme (26) will converge to the viscosity solution as the discretization is consistent, stable, and monotone. Hence, the following theorem follows from Lemmas 1-3.

Theorem 2.
e solution of the numerical scheme (26) converges to the viscosity solution of (10) as h ⟶ 0.

Iterative Algorithm for the Discrete System
In this section, to solve the nonlinear algebraic system (26) effectively, we propose and analyze an iterative algorithm for (26) at each time step, which is given in Algorithm 1. Now, we show the convergence of this iterative algorithm for the nonlinear algebraic system (26) in the following theorem.

Theorem 3. At each time step n of the iterative scheme, equation in Algorithm 1 generates a sequence of solutions V
l , starting from any initial guess V 0 , that converges monotonically to the solution V n+1 to (26) as l ⟶ ∞, provided that h is sufficiently small.
Proof. First, we show that the iterate V l is bounded for any l. Writing out the equation in Algorithm 1 in component form gives As in the proof of Lemma 2, it follows from β i ≥ 0, w i ≥ 0, and χ i ≤ 0 that Using a similar argument as used in proving Lemma 2, when h is sufficiently small, we have i.e., ‖V l ‖ ∞ is bounded independently of l.
Second, we show that the iterates V l form a nondecreasing sequence. e lth iteration of the equation in Algorithm 1 is given by which can be rewritten as follows: Subtracting (40) from the equation in Algorithm 1, we have e above scheme can be written as follows: Noting (42) us, Combining (42) with (45), we have (1) Let n � 0.
Further, from eorem 1, we know that I − (Δτ/2)S + ΔτP(V l ) is an M-matrix, so the inequality (46) implies by the discrete maximum principle, and hence, the iterates V l form a nondecreasing sequence.
To summarize, we have shown that the iterates V l form a nondecreasing sequence which is bounded from above. Hence, the sequence converges to the solution of numerical scheme (26). As for the uniqueness, suppose we have two solutions to the equation in Algorithm 1, B 1 and B 2 . en Some manipulation of (48a) and (48b) results in Using a similar argument as used in proving (46), we obtain Since I − (Δτ/2)S + ΔτP(B 2 ) is an M-matrix, it follows from (49) and (50) that By symmetry, we also have and hence B 1 � B 2 . is completes the proof.
Remark 4. From eorems 2 and 3, we know that the solution resulting from Algorithm 1 converges to the correct solution of (10), i.e., the unique viscosity solution or financially relevant solution.

Numerical Experiments
In this section, we present some numerical results to demonstrate the performance and convergence of the new numerical scheme. In the sequel, "error" refers to the difference between successive numerical solutions following mesh refinements, i.e., where V Δτ Δr denotes the computed solution on the spatial mesh Δr and time mesh size Δτ, and "order" is log 2 of the ratio of successive changes in option price as the grid is refined. Further, the penalty parameter c and the tolerance are chosen to be 10 3 and 10 − 8 in Algorithm 1, respectively. All codes were carried out in Matlab R2016a with 32.00 GB RAM and 3.20 GHz processor. e assumed parameters for American put options on a zero-coupon bond under the CIR model are (cf. [11]) as follows:  Tables 7-10 list the prices at θ � 0.08 and r � 0.9, 1.1, 1.3 for different values of σ, where r max � 3. From the results in these tables, we can observe that the accuracy of the numerical solutions of the American put options on a zero-coupon bond improves as the discretization grid is refined, and the order of convergence of the Crank-Nicolson scheme is roughly 2. is is consistent with the properties of the Crank-Nicolson scheme. Some authors have also used this approach when an analytic solution is not available [25,31,32]. Further, to visualize the numerical solutions, we plot the values of the American put options on a zero-coupon bond with different values of σ in Figures 1-3. From these figures, we can see that the numerical solutions are qualitatively excellent and the plots also indicate that our numerical scheme is stable.

Conclusion
We developed and analysed a Crank-Nicolson fitted finite volume method for pricing American options on a zerocoupon bond model. Based on a penalty approach, we first approximate the original problem by a nonlinear PDE. After that, a novel fitted finite volume method for the spatial discretization of the nonlinear penalized PDE system is coupled with the Crank-Nicolson time-stepping scheme. It is shown that this scheme is consistent, stable, and monotone, and hence, the convergence was guaranteed. To solve the nonlinear algebraic system effectively, an iterative algorithm was designed and its convergence was proved. Numerical experiments were performed to demonstrate the convergence, efficiency, and robustness of the new numerical method. Compared with the fully implicit time-stepping schemes proposed by Zhang et al. in [11,14], the Crank-Nicolson scheme here delivers more accurate approximate solutions.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest.