Pricing American bond options using a cubic spline collocation method

In this paper, American options on a discount bond are priced under the Cox-Ingrosll-Ross (CIR) model. The linear complementarity problem of the option value is solved numerically by a penalty method. The problem is transformed into a nonlinear partial differential equation (PDE) by adding a power penalty term. The solution of the penalized problem converges to the one of the original problem. To numerically solve this nonlinear PDE, we use the horizontal method of lines to discretize the temporal variable and the spatial variable by means of trapezoidal method and a cubic spline collocation method, respectively. We show that this full discretization scheme is second order convergent, and hence the convergence of the numerical solution to the viscosity solution of the continuous problem is guaranteed. Numerical results are presented and compared with other collocation methods given in the literature.


Introduction
Although the concept of interest rates seems to be something natural that everybody knows to deal with, the management of interest rate risk, i.e. the control of changes in future cash flows due to fluctuations in interest rates is an issue of great complexity.In particular, the pricing and hedging of products depending in large part on interest rates create the necessity for mathematical models.Here we adopt the CIR model developed by Cox, Ingersoll and Ross in 1985 [8], is one of the most widely used term structure model and has several favorable features.Various mathematical results about CIR were known since the work presented in [13].CIR model is defined as the unique strong solution of the equation 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 constant σ > 0. In practice, r is positive, which enforces the constraint 0 < σ 2 < 2κθ (Fellers condition [8]).In [8], it has been shown that the price P (r, t, s) of a pure discount bond with face value of one dollar at its maturity date s is given as follows P (r, t, s) = A(t, s)e −B(t,s)r , where and ζ is the market risk premium.We concentrate on the numerical solution of the CIR model which can be formulated as parabolic partial differential complementarity (PDC) problem with suitable boundary and terminal condition (see [1,17,22]).This complementarity problem is, in general, not analytically solvable.Hence, several approximation techniques have been developed for the solution of American bond option pricing problem.Numerical solution by penalty methods have been considered, e.g. by [1,2,9,12].In this paper we develop a novel numerical method for solving a PDC problem by using the cubic spline collocation method and the generalized Newton method.First, the PDC problem is approximated by a sequence of nonlinear equation problems by using the penalty method given in [23,27].It is shown that the penalized equation converges to the one of the PDC problem with an arbitrary order.This arbitrary order convergence rate allows one to achieve the required accuracy of the solution with a small penalty parameter.A numerical scheme for solving the penalized nonlinear PDE is also proposed.Then we apply the spline collocation method to approximate the solution of a boundary value problem of second order.The discret problem is formulated as to find the cubic spline coefficients of a nonsmooth system ϕ(Y ) = Y , where ϕ : R m → R m .In order to solve the nonsmooth equation we apply the generalized Newton method (see [4,5,20], for instance).We prove that the cubic spline collocation method converges quadratically provided that a property coupling the penalty parameter λ and the discretization parameter h is satisfied.
Pricing American bond options using a cubic spline collocation method 191 Numerical methods to approximate the solution of boundary value problems have been considered by several authors.We only mention the paper [3] and references therein, which use the spline collocation method for solving the boundary value problems.
The organization of this paper is as follows.In Section 2, we present the penalty method to approximate the PDC problem by a sequence of second order boundary value problems.In Section 3 we construct a cubic spline to approximate the solution of the boundary problem.Section 4 is devoted to the presentation of the generalized Newton method.In Section 5 we show the convergence of the cubic spline to the solution of the boundary problem and provide an error estimate.Finally, some numerical results are given in Section 6 to validate our methodology and compare our method with [25].

Penalty problem
Now, let u(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, the option pricing problem can be formulated as the following parabolic partial differential complementarity problem ( [24]): almost everywhere in (0, +∞) × (0, T ), where is a degenerate parabolic partial differential operator.The initial condition is u(r, 0) = Λ(r, 0) = max[K − P (r, T, s), 0] for a put, and the boundary conditions are For computational purposes, it is necessary to restrict r to a finite region where R denotes a sufficiently large number to ensure the accuracy of the solution, see [24].Thus, relation (2.2) becomes It is worth noting that T < s and K < P (0, T, s) = A(T, s) for a call option or K > A(T, s) for a put option, since otherwise the option would never be exercised and would be worthless.
Let λ be a real positive number.The penalty problem is given by the following boundary value problem (see [26]): where is the Black-Scholes operator defined as and with α(r) ≥ α > 0, γ(r) ≤ γ < 0 on Ω and α, β, γ, ψ 0 , ψ 1 , ψ 2 , are sufficiently smooth functions.Here we assume that the problem satisfies sufficient regularity and compatibility conditions which guarantee that the problem has a unique solution u ∈ C(Ω) C 2,1 (Ω) satisfying (see, [15,16,18]): where k is a constant.The following theorem proves the order of convergence of the solution u λ to u.
Theorem 2.2.( [14]) Let u and u λ be solutions to problems (2.1) and ( 2.3) respectively.Then, there exists a constant C > 0, independent of u , u λ and λ, such that where λ is the penalty parameter used in (2.3) Lemma 2.3.The function F is a nonlinear continuous and λ-Lipschitz on u λ .
Proof: F is Lipschitz, indeed For any two function u λ and v λ , we have: Pricing American bond options using a cubic spline collocation method 193 Then F is a nonlinear continuous function on u λ and satisfies the following Lipschitz condition:

✷
In the sequel of this paper we will consider the model (2.3) for numerical traitement.

Time discretization and description of the Trapezoidal method
We denote by .∞ the uniform norm.
The observation period has to be specified first and is set from today (time point t = 0) to time point T .This period will be divided into M equally spaced time intervals with length ∆t = T M , t m = m∆t.The price of the underlying (either the project value or the value of a twin security with the same risk profile) is assumed to stay in a range Ω r .We discretize the variable time in (2.3) by means of Trapezoidal method.Then the semi-discretization yields the following system of equations: That is where, L r : and u m λ approximate the exacte solution u λ (r, t) at the time level t m = m∆t.Then, the approximate problem of (2.3) is 2) where, for any m ≥ 0 and for any r ∈ Ω r , we have 2), at the (m + 1) th-time level.
The following theorem proves the order of convergence of the solution u m λ to u λ (r, t) Proof: We introduce the notation e m = u λ (r, t m ) − u m λ the error at step m.By Taylor series expansion of u λ , we have By using these expansions, we get Pricing American bond options using a cubic spline collocation method 195 and by Taylor series expansion of ∂u λ ∂t , we have )+O((∆t) 3 ), By using these expansions, and )), we have This implies By using this relation in (3.3) we get ))] +O((∆t) 3 ), by (3.1).Then, we obtain For an analytic function F we may bound the term O((∆t) 3 ) by c(∆t) 3 for some c > 0, and this upper bound is valid uniformly throughout [0, T ].Therefore, it follows from the Lipschitz condition of Lemma 2.3 and the triangle inequality that Clearly, the operator 1 ± ∆t 2 Lr satisfies a maximum principle (see, [6,7]) and consequently Since we are ultimately interested in letting ∆t → 0, there is no harm in assuming that ∆t.η < 2, with η = Lr ∞ + λ.We can thus deduce that We now claim that The proof is by induction on m.When m = 0 we need to prove that e0 ≤ 0 and hence that e0 = 0.This is certainly true, since at t0 = 0 the numerical solution matches the initial condition and the error is zero.
For general m ≥ 0, we assume that (3.5) is true up to m and use (3.4) to argue that This advances the inductive argument from m to m + 1 and proves that (3.5) is true.Since 0 < ∆t.η < 2, it is true that Consequently, relation (3.5) yields This bound is true for every nonnegative integer m such that m∆t < T .Therefore We deduce that u λ (r, tm) − u m λ ∞ ≤ C(∆t) 2 .In other words, problem (3.2) is second order convergent.✷ Pricing American bond options using a cubic spline collocation method 197 For any m ≥ 0, problem (3.2) has a unique solution and can be written on the following form: where It is easy to see that J λ is a nonlinear continuous function on u λ and satisfies the following Lipschitz condition: where In the sequel of this paper, we will focus on the solution of problem (3.6).

Spatial discretization and cubic spline collocation method
In this section we construct a cubic spline which approximates the solution u λ of problem (3.6), in the interval Ω r ⊂ R. We denote by .the Euclidean norm on R n+1 and S (k) the k th derivative of a function S. (Ω r , Θ) the space of piecewise polynomials of degree less than or equal to 3 over the subdivision Θ and of class C 2 everywhere on Ω r .Let B i , i = −3, • • • , n − 1, be the B-splines of degree 3 associated with Θ.These B-splines are positives and form a basis of the space S 4 (Ω r , Θ).
Consider the local linear operator Q 3 which maps the function u λ onto a cubic spline space S 4 (Ω r , Θ) and which has an optimal approximation order.This operator is the discrete C 2 cubic quasi-interpolant (see [19]) defined by where the coefficients µ j (u λ ) are determined by solving a linear system of equations given by the exactness of Q 3 on the space of cubic polynomial functions P 3 (Ω r ).Precisely, these coefficients are defined as follows: It is well known (see e.g.[11], chapter 5) that there exists constants C k , k = 0, 1, 2, 3, such that, for any function where u λ Ωr = max r∈Ωr |u λ (r)|.
By using the boundary conditions of problem (3.6), we obtain where From equation: (4.1), we can easily see that the spline S 1 satisfies the following equation The goal of this section is to compute a cubic spline collocation S λ = n−1 i=−3 c i,λ B i which satisfies the equation (3.6) at the points τ i , i = 0, ..., n + 2 with τ 0 = r 0 , and the coefficients c i,λ , i = −2, ..., n−2 satisfy the following collocation conditions: where ] T , and using equations (4.2) and (4.3), we get: Pricing American bond options using a cubic spline collocation method 199 and (P A (2) with where matrices A 0 , A 1 and A 2 are independent of h, with expression A 2 is: Results of this work are basically based on the invertibility of the matrix A 2 .Then, in order to prove that A 2 is invertible we give the flowing lemma.
j=−2 d j B j , then we have S(0) = S(R) = 0 and S (2) If we assume that S (2) = 0 in [r 0 , r n ], then using the Lemma 4.1 and the fact that S (2) has n + 1 zeros in [r 0 , r n ], we conclude that n + 1 ≤ n − 2, which is impossible.Therefore S (2) = 0 for each r ∈ Ω r .This means that the function S is a piecewise linear polynomial on Ω r .Since S(0) = S(R) = 0, then we obtain S(r) = 0 for any r ∈ Ω r .Consequently D = 0 and the matrix A 2 is invertible.✷ Then, relations (4.4) and (4.5) can be written in the following form with In order to determine the bounded of C − C ∞ , we need the following Lemma.
Proof: From the relation (4.8),We have For h sufficiently small, we conclude From the relation (4.9) and A 0 ∞ ≤ 1, we have For h sufficiently small, we conclude that h 2 ρ max i=0...n Pricing American bond options using a cubic spline collocation method 201 Proposition 4.4.Assume that the penalty parameter λ and the discretization parameter h satisfy the following relation: Then there exists a unique cubic spline which approximates the exact solution u λ of problem (3.6).
Proof: From (4.7), we have To prove the existence of cubic spline collocation it suffices to prove that ϕ admits a unique fixed point.Indeed, let Y 1 and Y 2 be two vectors of R n+1 .Then we have Using relation (3.7) and the fact that Then we obtain From relation (4.14), we conclude that Hence the function ϕ admits a unique fixed point.✷ In order to calculate the coefficients of the cubic spline collocation given by the nonsmooth system we propose the generalized Newton method defined by where I n+1 is the unit matrix of order n + 1 and V k is the generalized Jacobian of the function C λ → ϕ( C λ ), (see [4,5,20], for instance).

Convergence of the method
Theorem 5.1.Assume that the penalty parameter λ and the discretization parameter h satisfy the following relation: Then the cubic spline S λ converges to the solution u λ .Moreover the error estimate Proof: We pose ν = (I +U +V ) −1 (P A 2 ) −1 ∞ .From (4.8), (4.9) and Lemma 4.2, we have Since E is of order O( h 2 ∆t ), then there exists a constant K 1 such that Hence, we have On the other hand we have From relation (4.1), there exists a constant K 2 such that Using the fact that then, we obtain Pricing American bond options using a cubic spline collocation method 203 By using relation (5.2) and assumption (5.1) it is easy to see that (5.5) We have then from relations (5.3), (5.4) and (5.5), we deduce that Hence the proof is complete.✷ Remark 5.2.Theorem 5.1 provides a relation coupling the penalty parameter λ and the discretization parameter h, which guarantees the quadratic convergence of the cubic spline collocation S λ to the solution u λ of the penalty problem.

Numerical examples
In this section we verify numerically the obtained theoretical results in the previous section.If the exact solution is known, then at time t ≤ T the maximum error E max can be calculated as: Otherwise it can be estimated by the following double mesh principle: where S M,N (r, t) is the numerical solution on the M + 1 grids in space and N + 1 grids in time, and S 2M,2N (r, t) is the numerical solution on the 2M + 1 grids in space and 2N + 1 grids in time.Denote This last equation is nonlinear in C λ .We now apply the generalized Newton method to this equation.Note that when k > 1, we see that vector derivative of d at C λ satisfies d ′ ( C λ ) −→ ∞ as ψ − B C λ −→ 0 + .To overcome this difficulty, we use the technique proposed in [21] to smooth (6.1), yielding the following approximation to d( C λ ) : for k > 0, where 0 < ǫ << 1 is a transition parameter and W (z) is a function which smooths out the original d(z) around z = 0. We choose W (z) = a 1 + a 2 z + ... + a n z n−1 + a n+1 z n for n ≥ 3 and impose that W (z) is such that d(.) is smooth.This requires that W (z) satisties In this case, the function defined in (6.2) is globally smooth.Using the four conditions given in relation (6.1) and setting a 3 = ... = a n−1 = 0, we can easily find that Taking X = ψ − B C λ , then relation (6.2) can be expressed as In this section we apply the numerical method developed in this paper to an example presented in the paper of the authors Kai Zhang et al. [25].Then, we compare the obtained results to the ones given in [25].
A vanilla put option on a zero-coupon bond has the payoff Λ = max[K − P (r, T, s), 0], with P (r, T, s) = EA(T, s)e −B(T,s)r .
Pricing American bond options using a cubic spline collocation method 205 The parameters used for this put option on a bond under the CIR model are listed in Table 1.For the put option on a bond with the parameters in Table 1, we choose R = 2.The coarsest grid is defined as a uniform partition of the solution domain (0, 2) × (0, 1).
The comparison of the maximum error values between the method developed in this paper with the one developed in [25] will be taken at five different values of the number of space steps N = 201, 401, 801, 1601, 3201 and time steps M = 100, 200, 400, 800, 1600, for σ = 0.1 and σ = 0.5.We conduct experiments on different values of N , M and σ.Table 2 and Table 3 show values of the maximum error (max_error) obtained in our numerical experiments and the one obtained in [25].We see that the values of maximum error obtained by our method improve the ones in [25].

Conclusion
In this paper we have presented the American options on a bond under the Cox-Ingersoll-Ross model, this problem is approximated by a sequence of nonlinear equation problems by using the penalty method given in [23,27], and its time 0.00002 0.0024 discretization scheme.Then, we have developed and analyzed a cubic spline collocation method and the generalized Newton method for approximating solutions of the semi-discrete problem.We have shown the convergence of the method provided that the penalty and discrete parameters satisfy the relation (5.1).Moreover we have provided an error estimate of order O(h 2 ) with respect to the maximum norm .∞ .Numerical experiment was performed on one known model to validate the convergence and efficiency of the method.The computational results show that the proposed numerical method is an efficient alternative method to the one proposed in [25].

190A.
El hajaji, K. Hilal, A. Serghini and E. B. Mermri be a subdivision of the interval Ω r .Without loss of generality, we put r i = a + ih, where 0 ≤ i ≤ n and h = R n .Denote by S 4 (Ω r , Θ) := P 2 3

Table 1 :
Data used to value American put options on a zero-coupon bond under the CIR model.