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

Authors: Xiaoting Gan; Dengguo Xu
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

This paper develops and analyses a Crank–Nicolson fitted finite volume method 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 volume method 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.


---

## Body

## 1. 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–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–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. This 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–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 drift-dominated 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].The 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 first-order convergence rate. To the authors’ knowledge, based on the fitted finite volume method discretization, there is no theoretical result published about the Crank–Nicolson time-stepping 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. The numerical results further show that the Crank–Nicolson fitted finite volume scheme provides the rate of convergence is of approximately second order.The 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. Then, a Crank–Nicolson fitted finite volume scheme is constructed in Section3. 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. The paper ends with some conclusions.


## 2. 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:(1)dr=κθ−rdt+σrdW,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 σ2r is the variance with σ>0. In [24], it has been shown that the price Pr,t,s of a pure discount bond with face value $1 at its maturity date s is given as follows:(2)Pr,t,s=αt,se−βt,sr,where(3)αt,s=ϕ1eϕ2s−tϕ2eϕ1s−t−1+ϕ1ϕ3,βt,s=eϕ1s−t−1ϕ2eϕ1s−t−1+ϕ1,ϕ1=μ2+2σ2,ϕ2=μ+ϕ12,ϕ3=2θ¯σ2,θ¯=κθ,μ=κ+ζ,and ζ is the market risk premium. At the maturity date s, the price of a pure discount bond is its face value, i.e.,(4)Pr,t=s,s=1.Now, letVr,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 Vr,τ can be governed by the following PDCP (cf. [12]):(5)LVr,τ⩾0,Vr,τ−Λr,τ⩾0,LVr,τVr,τ−Λr,τ=0,a.e. in 0,+∞×0,T, where(6)LV=∂V∂τ−12σ2r∂2V∂r2−θ¯−μr∂V∂r+rV,and the initial condition is given by(7)Vr,0=Λr,0=maxK−Pr,T,s,0,and the boundary conditions are as follows:(8)V0,τ=Λr,τ,r⟶0,Vr,τ=Λr,τ,r⟶∞.For computational purpose, it is necessary to restrict the short-term interest rater in a finite region I=0,rmax, where rmax is a sufficiently large number to ensure the accuracy. Thus, the boundary conditions (8) become(9)V0,τ=Λ0,τ,Vrmax,τ=Λrmax,τ.By using the penalty method (cf. [11, 20, 25]), we can transform the PDCP (1) into the following nonlinear penalized PDE:(10)ℒVγS,τ−γΛ−Vγ+=0,with the given initial and boundary conditions:(11)Vγr,0=Λr,0,Vγ0,τ=K,VγR,τ=0,where γ is the penalty parameter satisfying γ>1 and x+=max0,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γ+ in (10) to be close to zero as γ becomes sufficiently large. Hence, the complementarity condition in (5) is approximately satisfied.Remark 1.
Since the degeneracy of the diffusion operatorsℒ, the nonsmoothness of payoff functions, and nonlinearity of max⋅,⋅, both problems (5) and (10) have no classical solutions, in general. Therefore, 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(12)Hx,u,Du,D2u=0inΩ,where H is a continuous function satisfying the ellipticity condition:(13)Hx,u,p,M≤Hx,u,p,NifM≥N,for any x∈Ω,Ω is a domain in RN, u∈R, p∈RN, if and only if ∀φ∈C2Ω, if x0∈Ω is a local maximum point of u−φ, one has(14)Hx0,ux0,Dφx0,D2φx0≤0,and ∀φ∈C2Ω, if x0∈Ω is a local minimum point of u−φ, one has(15)Hx0,ux0,Dφx0,D2φx0≥0.We note that, according to [11], the solution Vγ to penalized PDE (5) is an approximation to solution V of PDCP (1) in the viscosity sense. The convergence of Vγ to V as γ⟶∞ has been discussed by many researchers, e.g., Li and Wang [28] and Wang et al. [29]. Therefore, 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 γ of Vγ used in the above, but bear in mind that this V is the solution of penalized PDE (5).


## 3. 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).


### 3.1. Semidiscretization Scheme

Before proceeding to the Crank–Nicolson fitted finite volume discretization scheme, we first transform PDE (5) into the following conservative form:(16)∂V∂τ=∂∂rar∂V∂S+brV−crV+γΛ−V+,where(17)a=σ22>0,br=θ¯−μr−σ22,cr=r−μ.The fitted finite volume method is based on the self-adjoint form (16). Now we first define two spatial partitions of I=0,rmax. Let I be divided into M subintervals:(18)Ii=ri,ri+1,i=0,1,…,M−1,with 0=r0<r1<r2<⋅⋅⋅<rM=rmax. For each i=0,1,…,M−1, let hi=ri+1−ri. Also, we let ri−1/2=ri−1+ri/2andri+1/2=ri+ri+1/2 for each i=1,…,M−1. Then these intervals Ii∗=ri−1/2,ri+1/2, i=0,1,…,M, form a second partition of I=0,rmax if we define r−1/2=r0 and rM+1/2=rM.For eachi=1,…,M−1, integrating (16) over Ii∗, we have(19)∫Ii∗∂V∂τdr=ar∂V∂r+brVri−1/2ri+1/2−∫Ii∗crVdr+γ∫Ii∗Λ−V+dr.Applying the one-point quadrature rule to all the terms in (19) except the first term in the right-hand side, we obtain(20)∂Vi∂τli=ρVri+1/2−ρVri−1/2−ciliVi+γliΛi−Vi+,i=1,…,M−1,where li=ri+1/2−ri−1/2 is the length of interval Ii∗ and ρV is the weighted flux density associated with V defined as follows:(21)ρV≔arV′+brV,where Vi and Λi denote the approximation of Vri,τ and Λri,τ, 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(22)∂V∂τ=SV+PΛ−V+F,where V=V1,V2,…,VM−1T, Λ=Λ1,Λ2,…,ΛM−1T, and P is a M−1×M−1 diagonal matrix with(23)Pii=γ,ifVi<Λi,0,otherwise,and S is a M−1×M−1 semidiscretization matrix given by(24)S=χ1w1β2χ2w2⋱⋱⋱βM−1χM−1M−1×M−1,where(25)β1=12l1a−b1/2,w1=b3/2r2η1r2η1−r1η1l1,χ1=−12l1a+b1/2−b3/2r1η1r2η1−r1η1l1−c1,βi=bi−1/2ri−1ηi−1riηi−1−ri−1ηi−1li,χi=−bi−1/2riηi−1riηi−1−ri−1ηi−1li−bi+1/2riηiri+1ηi−riηili−ci,wi=bi+1/2ri+1ηiri+1ηi−riηili,for i=2,…,M−1, F=β1V0,0,…,0,wM−1VMT, and V0 and VM being equal to the given boundary conditions in (9), and ηi=bi+1/2/a.


### 3.2. The Crank–Nicolson Scheme

We now consider the time-stepping 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. Thus, the Crank–Nicolson scheme of (22) can be written in the following matrix form:(26)I−Δτ2S+ΔτPVn+1Vn+1=I+Δτ2SVn+ΔτPVn+1Λ+ΔτF,n.for n=0,1,…,N−1, where Vn=V1n,V2n,…,VM−1nT, Fn=β1V0n,0,…,0,wM−1VMnT, and PVn+1 is a M−1×M−1 diagonal matrix with(27)PVn+1ii=γ,ifVin+1<Λi,0,otherwise,where I is a M−1×M−1 unit matrix and Vin denotes the approximations of Vri,τ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 PVn+1, we further have the following result for the numerical scheme (26).Theorem 1.
Both matricesI−Δτ/2S and I−Δτ/2S+ΔτPVn+1 are M-matrices.Remark 2.
A realn×n matrix A=aij with aij≤0 for all i≠j is an M-matrix if A is nonsingular and A−1≥O, and further if A=aij is a real, irreducibly diagonally dominant n×n matrix with aij≤0 for all i≠j, and aii>0 for all 1≤i≤n, then A−1≥O (cf. Varga [30], p. 91).Remark 3.
Theorem1 implies that the Crank–Nicolson scheme (26) satisfies the discrete maximum principle, which guarantees that the discrete arbitrage inequality holds in option pricing theory.


## 3.1. Semidiscretization Scheme

Before proceeding to the Crank–Nicolson fitted finite volume discretization scheme, we first transform PDE (5) into the following conservative form:(16)∂V∂τ=∂∂rar∂V∂S+brV−crV+γΛ−V+,where(17)a=σ22>0,br=θ¯−μr−σ22,cr=r−μ.The fitted finite volume method is based on the self-adjoint form (16). Now we first define two spatial partitions of I=0,rmax. Let I be divided into M subintervals:(18)Ii=ri,ri+1,i=0,1,…,M−1,with 0=r0<r1<r2<⋅⋅⋅<rM=rmax. For each i=0,1,…,M−1, let hi=ri+1−ri. Also, we let ri−1/2=ri−1+ri/2andri+1/2=ri+ri+1/2 for each i=1,…,M−1. Then these intervals Ii∗=ri−1/2,ri+1/2, i=0,1,…,M, form a second partition of I=0,rmax if we define r−1/2=r0 and rM+1/2=rM.For eachi=1,…,M−1, integrating (16) over Ii∗, we have(19)∫Ii∗∂V∂τdr=ar∂V∂r+brVri−1/2ri+1/2−∫Ii∗crVdr+γ∫Ii∗Λ−V+dr.Applying the one-point quadrature rule to all the terms in (19) except the first term in the right-hand side, we obtain(20)∂Vi∂τli=ρVri+1/2−ρVri−1/2−ciliVi+γliΛi−Vi+,i=1,…,M−1,where li=ri+1/2−ri−1/2 is the length of interval Ii∗ and ρV is the weighted flux density associated with V defined as follows:(21)ρV≔arV′+brV,where Vi and Λi denote the approximation of Vri,τ and Λri,τ, 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(22)∂V∂τ=SV+PΛ−V+F,where V=V1,V2,…,VM−1T, Λ=Λ1,Λ2,…,ΛM−1T, and P is a M−1×M−1 diagonal matrix with(23)Pii=γ,ifVi<Λi,0,otherwise,and S is a M−1×M−1 semidiscretization matrix given by(24)S=χ1w1β2χ2w2⋱⋱⋱βM−1χM−1M−1×M−1,where(25)β1=12l1a−b1/2,w1=b3/2r2η1r2η1−r1η1l1,χ1=−12l1a+b1/2−b3/2r1η1r2η1−r1η1l1−c1,βi=bi−1/2ri−1ηi−1riηi−1−ri−1ηi−1li,χi=−bi−1/2riηi−1riηi−1−ri−1ηi−1li−bi+1/2riηiri+1ηi−riηili−ci,wi=bi+1/2ri+1ηiri+1ηi−riηili,for i=2,…,M−1, F=β1V0,0,…,0,wM−1VMT, and V0 and VM being equal to the given boundary conditions in (9), and ηi=bi+1/2/a.


## 3.2. The Crank–Nicolson Scheme

We now consider the time-stepping 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. Thus, the Crank–Nicolson scheme of (22) can be written in the following matrix form:(26)I−Δτ2S+ΔτPVn+1Vn+1=I+Δτ2SVn+ΔτPVn+1Λ+ΔτF,n.for n=0,1,…,N−1, where Vn=V1n,V2n,…,VM−1nT, Fn=β1V0n,0,…,0,wM−1VMnT, and PVn+1 is a M−1×M−1 diagonal matrix with(27)PVn+1ii=γ,ifVin+1<Λi,0,otherwise,where I is a M−1×M−1 unit matrix and Vin denotes the approximations of Vri,τ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 PVn+1, we further have the following result for the numerical scheme (26).Theorem 1.
Both matricesI−Δτ/2S and I−Δτ/2S+ΔτPVn+1 are M-matrices.Remark 2.
A realn×n matrix A=aij with aij≤0 for all i≠j is an M-matrix if A is nonsingular and A−1≥O, and further if A=aij is a real, irreducibly diagonally dominant n×n matrix with aij≤0 for all i≠j, and aii>0 for all 1≤i≤n, then A−1≥O (cf. Varga [30], p. 91).Remark 3.
Theorem1 implies that the Crank–Nicolson scheme (26) satisfies the discrete maximum principle, which guarantees that the discrete arbitrage inequality holds in option pricing theory.


## 4. 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.Lemma 1.
The numerical 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 uh be the discrete approximation of u. It has been shown in [11, 16] that the exact and the discrete flux yield(28)ρu−ρhuhri+1/2=Oh,where ρhV satisfies(29)ρhV=ρiV,∀r∈Ii,for i=0,1,…,M−1, and h=maxmaxihi,Δτ be the mesh parameter. Hence, it is easy to see that the consistency of the numerical scheme (26) is the consequent result.The stability result for the numerical scheme (11) is given as follows.Lemma 2.
The numerical scheme (26) is stable as h⟶0, i.e.,(30)Vn∞≤Λ∞.Proof.
Writing out (26) in component form gives(31)1−Δτ2χi+ΔτPVn+1iiVin+1−Δτ2βiVi−1n+1−Δτ2wiVi+1n+1=1+Δτ2χiVin+Δτ2βiVi−1n+Δτ2wiVi+1n+ΔτPVn+1iiΛi+ΔτFin.
It follows fromβi≥0,wi≥0,χi≤0 that, when h is sufficiently small(32)1−Δτ2χi+ΔτPVn+1iiVin+1≤Δτ2βiVi−1n+1+Δτ2wiVi+1n+1+1+Δτ2χiVin+Δτ2βiVi−1n+Δτ2wiVi+1n+ΔτPVn+1iiΛi+ΔτFin≤Δτ2βi+wiVn+1∞+1+Δτ2χi+βi+wiVn∞+ΔτPVn+1iiΛ∞+ΔτFn∞.
Choosel∈1,…,M−1 such that Vn+1∞=Vln+1. Then, (32) becomes(33)1−Δτ2χl+ΔτPVn+1llVn+1∞≤Δτ2βl+wlVn+1∞+1+Δτ2χl+βl+wlVn∞+ΔτPVn+1llΛ∞.
Then, we obtain the following estimate:(34)Vn+1∞≤maxVn∞,Λ∞×1+Δτ/2χl+βl+wl+ΔτPVn+1ll1+Δτ/2−χl−βl−wl+ΔτPVn+1ll=maxVn∞,Λ∞×1−Δτ/2rl+ΔτPVn+1ll1+Δτ/2rl+ΔτPVn+1ll≤maxVn∞,Λ∞,where −χl−βl−wl=rl>0. Combining (7) and (34) gives(35)Vn+1∞≤maxVn∞,Λ∞≤⋯≤maxV0∞,Λ∞=Λ∞,for all feasible n. Hence, the numerical scheme (26) is stable.The monotonicity of the numerical scheme (26) is given in the following lemma.Lemma 3.
The numerical scheme (26) is unconditionally monotone.Proof.
Fori=0 or M, the lemma is trivially true. When 0<i<M, it follows from Theorem 1 that the matrices −Δτ/2S and I−Δτ/2S are M-matrices, and hence, I−Δτ/2SVn+1i is a strictly increasing function of Vin+1 and nonincreasing function of Vi+1n+1 and Vi−1n+1, and −Δτ/2SVni is also a strictly increasing function of Vin and nonincreasing function of Vi+1n and Vi−1n. On the other hand, it is obvious that −Vin is a decreasing function of Vin. Finally, the penalty term −ΔτPVn+1iiΛi−Vin+1 is a nondecreasing function of Vin+1. 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.
The solution of the numerical scheme (26) converges to the viscosity solution of (10) as h⟶0.


## 5. 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.Algorithm 1: Iterative algorithm for (26).
(1)
Letn=0.(2)
Setl=0 and V^0=Vn.(3)
SolveI−Δτ/2S+ΔτPV^lV^l+1=I+Δτ/2SVn+ΔτPV^lΛ+ΔτFn(4)
IfmaxiV^il+1−V^il/max1,V^il+1<tolerance then stop; otherwise, set l:=l+1 and go to Step 3.(5)SetVn+1=V^l,n=n+1 and go to Step 2.Now, we show the convergence of this iterative algorithm for the nonlinear algebraic system (26) in the following theorem.Theorem 3.
At each time stepn 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 Vn+1 to (26) as l⟶∞, provided that h is sufficiently small.Proof.
First, we show that the iterateV^l is bounded for any l. Writing out the equation in Algorithm 1 in component form gives(36)1−Δτ2χi+ΔτPV^liiV^il+1−Δτ2βiV^i−1l+1−Δτ2wiV^i+1l+1=1+Δτ2χiVin+Δτ2βiVi−1n+Δτ2wiVi+1n+ΔτPV^liiΛi+ΔτFin.
As in the proof of Lemma2, it follows from βi≥0,wi≥0, and χi≤0 that(37)1−Δτ2χi+ΔτPqV^liiV^il+1≤Δτ2βiVi−1l+1+Δτ2wiVi+1l+1+1+Δτ2χiVin+Δτ2βiVi−1n+Δτ2wiVi+1n+ΔτPVliiΛi+ΔτFin≤Δτ2βi+wiV^l+1∞+1+Δτ2χi+βi+wiVn∞+ΔτPV^liiΛ∞+ΔτFn∞.
Using a similar argument as used in proving Lemma2, when h is sufficiently small, we have(38)V^l+1∞≤Λ∞,i.e., V^l∞ is bounded independently of l.
Second, we show that the iteratesV^l form a nondecreasing sequence. The lth iteration of the equation in Algorithm 1 is given by(39)I−Δτ2S+ΔτPV^l−1V^l=I+Δτ2SVn+ΔτPV^l−1Λ+ΔτFn,which can be rewritten as follows:(40)I−Δτ2SV^l=I+Δτ2SVn+ΔτPV^l−1Λ−V^l+ΔτFn.
Subtracting (40) from the equation in Algorithm 1, we have(41)I−Δτ2SV^l+1−V^l=ΔτPV^l−PV^l−1Λ+ΔτPV^l−1V^l−PV^lV^l+1.
The above scheme can be written as follows:(42)I−Δτ2S+ΔτPV^lV^l+1−V^l=ΔτPV^l−PV^l−1Λ−V^l.
Noting (42), we examine each component of PV^l−PV^l−1Λ−V^l in the following two cases.Case 1.
IfV^il<Λi, then we have PV^lii=γ, and hence(43)PV^l−PV^l−1Λ−V^l=γ−PV^l−1iiΛi−V^il≥0.Case 2.
IfV^il≥Λi, then we have PV^lii=0, and hence(44)PV^l−PV^l−1Λ−V^l=−PV^l−1iiΛi−V^il≥0.
Thus,(45)PV^l−PV^l−1Λ−V^l≥0.
Combining (42) with (45), we have(46)I−Δτ2S+ΔτPV^lV^l+1−V^l≥0.
Further, from Theorem1, we know that I−Δτ/2S+ΔτPV^l is an M-matrix, so the inequality (46) implies(47)V^l+1−V^l≥0,by the discrete maximum principle, and hence, the iterates V^l form a nondecreasing sequence.
To summarize, we have shown that the iteratesV^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 Algorithm1, B1 and B2. Then(48a)I−Δτ2S+ΔτPB1B1=I+Δτ2SVn+ΔτPB1Λ+ΔτFn,(48b)I−Δτ2S+ΔτPB2B2=I+Δτ2SVn+ΔτPB2Λ+ΔτFn.
Some manipulation of (48a) and (48b) results in(49)I−Δτ2S+ΔτPB2B1−B2=ΔτPB1−PB2Λ−B1.
Using a similar argument as used in proving (46), we obtain(50)PB1−PB2Λ−B1≥0.
SinceI−Δτ/2S+ΔτPB2 is an M-matrix, it follows from (49) and (50) that(51)B1≥B2.
By symmetry, we also have(52)B2≥B1,and hence B1=B2. This completes the proof.Remark 4.
From Theorems2 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.


## 6. 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.,(53)error=VΔrΔτ−VΔr/2Δτ/2,where VΔrΔτ denotes the computed solution on the spatial mesh Δr and time mesh size Δτ, and “order” is log2 of the ratio of successive changes in option price as the grid is refined. Further, the penalty parameter γ and the tolerance are chosen to be 103 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.The assumed parameters for American put options on a zero-coupon bond under the CIR model are (cf. [11]) as follows:(54)κ=0.1,θ=0.03or0.08,σ=0.1,0.3,0.5or0.8,ζ=0,E=100,K=60,T=1,s=5,rmax=2or3,and then the grid is defined as consecutive uniform partitions of the domain 0,rmax×0,T with m and n being the space and time steps, respectively.Tables1-2 list the prices of the American put option on a zero-coupon bond under the CIR model at θ=0.03 and r=0.9,1.0,1.1 for different values of σ, where rmax=2. Tables 3–6 list the prices at θ=0.08 and r=0.8,1.0,1.5 for different values of σ, where rmax=2. And Tables 7–10 list the prices at θ=0.08 and r=0.9,1.1,1.3 for different values of σ, where rmax=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. This 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.Table 1
The prices of the American put option on a zero-coupon bond withθ=0.03 and σ=0.1 as m=n, where rmax=2.
mr=0.9r=1.0r=1.1ValueErrorOrderValueErrorOrderValueErrorOrder10054.581754814——56.055386819——57.120900432——20054.5817569922.18e − 6—56.0553886701.85e − 6—57.1209019711.54e − 6—40054.5817575385.46e − 72.056.0553891344.64e − 72.057.1209023563.85e − 72.080054.5817576741.36e − 72.056.0553892501.16e − 72.057.1209024529.64e − 82.0160054.5817577093.41e − 82.056.0553892792.90e − 82.057.1209024762.41e − 82.0Table 2
The prices of the American put option on a zero-coupon bond withθ=0.03 and σ=0.8 as m=n, where rmax=2.
mr=0.9r=1.0r=1.1ValueErrorOrderValueErrorOrderValueErrorOrder10036.304526470——39.768020936——42.722357080——20036.3045278031.33e − 6—39.7680222231.29e − 6—42.7223583041.22e − 6—40036.3045281363.33e − 72.039.7680225443.22e − 72.042.7223586103.06e − 72.080036.3045282198.33e − 82.039.7680226248.04e − 82.042.7223586867.65e − 82.0160036.3045282402.08e − 82.039.7680226452.01e − 82.042.7223587061.91e − 82.0Table 3
The prices of the American put option on a zero-coupon bond withθ=0.08 and σ=0.1 as m=n, where rmax=2.
mr=0.8r=1.0r=1.5ValueErrorOrderValueErrorOrderValueErrorOrder10052.797853729——56.188221908——59.162911163——20052.7978557021.97e − 6—56.1882234721.56e − 6—59.1629117395.76e − 7—40052.7978561974.95e − 72.056.1882238633.91e − 72.059.1629118831.44e − 72.080052.7978563201.23e − 72.056.1882239619.80e − 82.059.1629119193.60e − 82.0160052.7978563513.10e − 82.056.1882239852.40e − 82.059.1629119289.00e − 92.0Table 4
The prices of the American put option on a zero-coupon bond withθ=0.08 and σ=0.3 as m=n, where rmax=2.
mr=0.8r=1.0r=1.5ValueErrorOrderValueErrorOrderValueErrorOrder10049.640003215——54.020836795——58.432510559——20049.6400044321.22e − 6—54.0208377579.62e − 7—58.4325109624.03e − 7—40049.6400047363.04e − 72.054.0208379982.41e − 72.058.4325110621.00e − 72.080049.6400048127.60e − 82.054.0208380586.00e − 82.058.4325110872.50e − 82.0160049.6400048311.90e − 82.054.0208380731.50e − 82.058.4325110936.00e − 92.0Table 5
The prices of the American put option on a zero-coupon bond withθ=0.08 and σ=0.5 as m=n, where rmax=2.
mr=0.8r=1.0r=1.5ValueErrorOrderValueErrorOrderValueErrorOrder10043.763440286——49.554320053——56.486342717——20043.7634416741.39e − 6—49.5543212451.19e − 6—56.4863433546.37e − 7—40043.7634420213.47e − 72.049.5543215432.98e − 72.056.4863435131.59e − 72.080043.7634421088.70e − 82.049.5543216177.40e − 82.056.4863435534.00e − 82.0160043.7634421302.20e − 82.049.5543216361.90e − 82.056.4863435631.00e − 82.0Table 6
The prices of the American put option on a zero-coupon bond withθ=0.08 and σ=0.8 as m=n, where rmax=2.
mr=0.8r=1.0r=1.5ValueErrorOrderValueErrorOrderValueErrorOrder10032.874384467——40.226739920——50.994927357——20032.8743856791.21e − 6—40.2267411341.21e − 6—50.9949282358.78e − 7—40032.8743859843.05e − 72.040.2267414373.03e − 72.050.9949284552.20e − 72.080032.8743860617.70e − 82.040.2267415137.60e − 82.050.9949285105.50e − 82.0160032.8743860801.90e − 82.040.2267415321.90e − 82.050.9949285241.40e − 82.0Table 7
The prices of the American put option on a zero-coupon bond withθ=0.08 and σ=0.1 as m=n, where rmax=3.
mr=0.9r=1.1r=1.3ValueErrorOrderValueErrorOrderValueErrorOrder15054.765185055——57.217095894——58.497251824——30054.7651868491.79e − 6—57.2170972181.32e − 6—58.4972527198.95e − 7—60054.7651872994.49e − 72.057.2170975493.31e − 72.058.4972529432.24e − 72.0120054.7651874111.12e − 72.057.2170976328.29e − 82.058.4972529995.61e − 82.0240054.7651874392.81e − 82.057.2170976532.07e − 82.058.4972530131.40e − 92.0Table 8
The prices of the American put option on a zero-coupon bond withθ=0.08 and σ=0.3 as m=n, where rmax=3.
mr=0.9r=1.1r=1.3ValueErrorOrderValueErrorOrderValueErrorOrder15052.133287029——55.449461279——57.347886873——30052.1332881261.10e − 6—55.4494621078.28e − 7—57.3478874615.88e − 7—60052.1332884002.74e − 72.055.4494623142.07e − 72.057.3478876081.47e − 72.0120052.1332884686.85e − 82.055.4494623655.17e − 82.057.3478876453.68e − 82.0240052.1332884851.71e − 82.055.4494623781.29e − 82.057.3478876549.19e − 92.0Table 9
The prices of the American put option on a zero-coupon bond withθ=0.08 and σ=0.5 as m=n, where rmax=3.
mr=0.9r=1.1r=1.3ValueErrorOrderValueErrorOrderValueErrorOrder15046.979531081——51.615448216——54.585318188——30046.9795323821.30e − 6—51.6154492901.07e − 6—-54.5853190308.41e − 7—60046.9795327073.25e − 72.051.6154495582.68e − 72.054.5853192402.10e − 72.0120046.9795327888.13e − 82.051.6154496266.71e − 82.054.5853192925.26e − 82.0240046.9795328092.03e − 82.051.6154496421.68e − 82.054.5853193061.31e − 82.0Table 10
The prices of the American put option on a zero-coupon bond withθ=0.08 and σ=0.8 as m=n, where rmax=3.
mr=0.9r=1.1r=1.3ValueErrorOrderValueErrorOrderValueErrorOrder15036.842142692——43.113757292——47.676423742——30036.8421439401.25e − 6—43.1137584531.16e − 6—47.6764247681.03e − 6—60036.8421442523.12e − 72.043.1137587432.90e − 72.047.6764250252.57e − 72.0120036.8421443307.80e − 82.043.1137588167.26e − 82.047.6764250896.42e − 82.0240036.8421443491.95e − 82.043.1137588341.81e − 82.047.6764251051.60e − 82.0Figure 1
The American bond option value surface and the option value atτ=T based on m,n=400,400, where rmax=2. (a) θ = 0.03, σ = 0.1; (b) θ = 0.03, σ = 0.8; (c) θ = 0.03, σ = 0.1; (d) θ = 0.03, σ = 0.8.
(a)(b)(c)(d)Figure 2
The American bond option value surface and the option value atτ=T based on m,n=400,400, where rmax=2. (a) θ = 0.08, σ = 0.1; (b) θ = 0.08, σ = 0.3; (c) θ = 0.08, σ = 0.1; (d) θ = 0.08, σ = 0.3.
(a)(b)(c)(d)Figure 3
The American bond option value surface and the option value atτ=T based on m,n=600,600, where rmax=3. (a) θ = 0.08, σ = 0.5; (b) θ = 0.08, σ = 0.8; (c) θ = 0.08, σ = 0.5; (d) θ = 0.08, σ = 0.8.
(a)(b)(c)(d)


## 7. Conclusion

We developed and analysed a Crank–Nicolson fitted finite volume method for pricing American options on a zero-coupon 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.



---

*Source: 1052084-2020-05-28.xml*