Comparison of numerical methods to price zero coupon bonds in a two-factor CIR model

In this paper we price a zero coupon bond under a Cox–Ingersoll–Ross (CIR) two-factor model using various numerical schemes. To the best of our knowledge, a closed-form or explicit price functional is not trivial and has been less studied. The use and comparison of several numerical methods to determine the bond price is one contribution of this paper. Ordinary differential equations (ODEs) , finite difference schemes and simulation are the three classes of numerical methods considered. These are compared on the basis of computational efficiency and accuracy, with the second aim of this paper being to identify the most efficient numerical method. The numerical ODE methods used to solve the system of ODEs arising as a result of the affine structure of the CIR model are more accurate and efficient than the other classes of methods considered, with the Runge–Kutta ODE method being the most efficient. The Alternating Direction Implicit (ADI) method is the most efficient of the finite difference scheme methods considered, while the simulation methods are shown to be inefficient. Our choice of considering these methods instead of the other known and apparently new numerical methods (eg Fast Fourier Transform (FFT) method, Cosine (COS) method, etc.) is motivated by their popularity in handling interest rate instruments.


INTRODUCTION
1.1 Equilibrium term-structure models like the ones we consider in this paper (two-factor Cox-Ingersoll-Ross (CIR) models) are common in non-defaultable bond pricing. Since a coupon-paying bond can be stripped into sums of zero-coupon bonds, the focus has always been on pricing zero-coupon bonds. With affine term structures, the expectation is to get an affine yield curve whose components can be explicitly obtained. There are other cases where closed-form pricing solutions of such bonds may not be either easy to obtain (like in our case) or almost impossible to get using classical methods. Our paper will not concentrate on justifying the right term-structure model, but on comparison of numerical methods used to price the zero coupon bond given a CIR two-factor model structure. To the best of our knowledge this is the first time such rigorous numerical comparison and subsequent simulation using real-life data has been done.

1.2
The model considered in this paper is a time-homogeneous short-rate equilibrium two-factor model. The two most prominent short-rate models are arguably the single-factor models proposed by Vasicek (1977) and Cox et al. (1985). Arguments for more than one factor models can be found, for example, in Stambaugh (1988) and Chen & Scott (2003) who argue that two-or three-factor models are required to adequately represent market dynamics. Litterman & Scheinkman (1991) propose the use of three-factor models while Jamshidian & Zhu (1996) argue that a single factor explains approximately 70 per cent of the variation in yield curves, with a second factor explaining approximately 15 per cent. The focus of this paper is on pricing a zero-coupon bond in a two-factor extension of the CIR model by means of numerical methods.

1.3
As mentioned, to the best of our knowledge, a closed-form solution to the price of the non-defaultable zero-coupon bond in this term-structure setup has not been dealt with adequately. We fill the gap by showing that several numerical methods adequately give best approximations and this is shown in the numerical results provided. The use of numerical methods by Duffie & Kan (1996) to estimate the solution of a partial differential equation (PDE) of a similar form to that presented in this paper and the argument by Brennan & Schwartz (1979) that no closed-form solution exists for a similar PDE is reason to believe that a closed-form solution does not exist. This allows for the identification of efficient methods that can be used to construct yield curves.

1.4
Literature comparing some of the numerical methods used in this paper is considered in section 3, although few comparisons of applications similar to those presented in this paper were found. One of the contributions of this paper is the comparison of several numerical methods across three different classes of methods: ordinary differential equation (ODE) methods used to solve the system of ODEs arising as a result of the affine structure of the CIR model; finite difference scheme methods applied to the PDE of the CIR model; and simulation methods applied to the equation of the bond price. A second contribution of this paper is the identification of the most efficient numerical method of all those considered. The results in this paper may also be relevant to other short-rate models, although not all classes of numerical methods are applicable to other short-rate models.

1.5
This paper begins with a description of the two-factor CIR model in section 2. This section also contains the derivation of three equations, any of which can be solved to obtain the price of zero-coupon bonds in the two-factor CIR model. Section 3 describes the numerical methods considered in this paper, after which a description of the methodology used for the comparison is described in section 4. Section 5 contains the results of the comparison, with section 6 containing the conclusions.

2.1
Without loss of generality, the two-factor CIR model considered in this paper is the specification of the model according to Shreve (2004), and is given by the short rate under the risk neutral measure Q, where . By construction, we expect that in probability, ( ) 0 R t ≥ for all 0 t ≥ , ensuring consistency with the market dynamic of positive interest rates. The above CIR model is in canonical form (specified such that it contains the least number of parameters), from which more complex two-factor affine yield models can be derived (Shreve, 2004).

2.2
Throughout this paper it is assumed the non-defaultable zero-coupon bond (simply referred to as bond in this paper) will mature at time T with redemption of 1, with the bond price being calculated at the current time t T ≤ . The remaining time to maturity is T t τ = − . The bond price is denoted as ( ) , P t T and since the price process is Markovian, we have ( ) ( ) 1 2 , , , P t T f t y y = for some measurable function f.

2.3
In the rest of the paper, we shall use the notation f x to imply f x ∂ ∂ . Then the Black-Scholes partial differential equation (PDE)   , , 1 f T y y = for all y 1 ≥ 0, y 2 ≥ 0 and for all 0 t T ≤ < . The above PDE is linear and of order two in three variables. The bond price ( ) Y t y = ) can be calculated directly (where possible) by finding closedform solutions to this PDE given the unique boundary condition. Where that is not easy or possible, like in this case, we use numerical solutions to PDEs applied to eq. (2). We shall be using finite difference methods as part of solution methods in the next section.

2.4
Pricing this bond could be done in many different ways. One way is through simulation (eg Monte Carlo simulation) of eq. (1). As stated before, if a closed-form solution to the price exists, we could, as another method, explicitly solve the PDE eq. (2). That could easily be achieved, seeing that the bond price can be specified as: for some functions ( ) ( )

2.5
Solving for ( ) ( ) The initial conditions are 2.6 We posit that the non-linear system of ODEs given by eq. (4) are not trivial in finding explicit solutions. We are less interested in the resulting phase planes and thus we in part use the ODE class of numerical methods to estimate the bond price apart from the finite difference methods proposed before. We show that our choice of numerical methods works perfectly because all give good approximations to at least three decimal places.

3.1
Simulation methods, finite difference scheme (FD) methods and ODE methods are the three categories of numerical methods considered in this paper. Multiple methods within each category are considered. The two simulation methods considered are Monte Carlo simulation and antithetic variates. The five FD methods considered are the alternating direction implicit (ADI) method, the explicit method, the Hopscotch method, the implicit method and the Crank-Nicolson (CN) method. The five numerical methods considered for ODE are the implicit method, the Euler/explicit method, the Runge-Kutta method, the Taylor method (second order) and the Crank-Nicolson (CN) method. Note that the last group target the nonlinear system of ODEs (4) while the FD methods are used for the PDE boundary value eq. (2). The simulation methods are used in eq. (1). We show by way of results for chosen parameters how the three sets of methods give approximate values for ( ) 1 2 , , f t y y .

3.2
Little prior research has been found on the comparison of numerical methods in pricing non-defaultable zero-coupon bonds across the three classes of methods considered, with this comparison being one of the contributions of this paper. An immediate comparison between FD methods and simulation methods can be found in Boyle (1977) who proposes that simulation methods are computationally inefficient but have the advantage of flexibility while Wilmott (2006) argues that for models containing less than four random factors, finite difference schemes are more efficient than simulation methods. We shall look more closely into these comparisons in section 5.

3.3
Few examples have been found in prior research of comparisons of some of the FD methods considered in this paper in an application that is similar to the one in this paper. Geske & Shastri (1985) compare the computational efficiency of the implicit and explicit methods when used to solve the Black-Scholes PDE, showing the explicit method to take roughly 60 per cent of the computation time taken by the implicit method. Hull & White (1990) argue that the explicit scheme uses between 40 and 70 per cent of the computation time of the implicit method. Although not in an application similar to that in this paper, previous results when using these methods (see eg Cairns, 2004) show that the CN method converges more quickly than the explicit and implicit methods, while the rate of convergence of the implicit method may be quicker than the explicit method in certain applications, but not all. Duffie & Kan (1996) argues that the ADI method is less computationally intensive than the standard implicit method and the Crank-Nicolson method, but is not guaranteed to result in stable solutions. They further showed that the errors resulting from the ADI method are small when used to solve a three dimensional PDE similar to that in this paper. The results from some of theses papers are compared to the results in this paper in section 5.

3.4
Simulation methods 3.4.1 We briefly summarise the Monte Carlo simulation. For more detailed treatment, the reader can refer, for example, to Cairns (2004:185) and references therein.

The time interval
and 0 k n ≤ ≤ . A simulated path of the short rate is obtained as where 1k z and 2k z are independent random simulations of the standard normal distribution.
is the estimated bond price.

3.4.3
The antithetic variates method is a variance reduction method using a similar procedure to that specified above. Each iteration involves two simulated paths ( ) are obtained, labelled A i1 and A i2 . The first simulated path uses the values z 1i and z 2i as specified above, the second uses -z 1i and -z 2i . The estimate A i in each iteration is calculated as the average of A i1 and A i2 . 3.4.4 There are two causes of error in the estimated price ( ) , P t T . The first is the variation in the sample estimate of the expectation, reduced as the number of simulations increases. The second is the discretisation error of the short rate and the error in estimating an integral with a discrete sum, reduced as Δ t →0. The antithetic variates method is expected to result in a lower variance of the estimated price than Monte Carlo simulation (Cairns, 2004). Boyle et al. (1997), however, argues that the antithetic variates method is inefficient in reducing the variance of the estimated price, and suggests the alternatives of control variates, Latin hyper cube sampling and moment matching. Joy et al. (1996) suggests quasi-Monte Carlo methods as an additional alternative. These alternative methods are not considered in this paper.

3.5
Finite difference scheme (FD) methods 3.5.1 ovErviEW 3.5.1.1 In order to specify eq. (2) in terms of bounded variables for easier application of the FD methods, y 1 and y 2 were transformed to the variables θ 1 and θ 2 . The bond price is subsequently represented as ( ) f t θ θ . The common transformation used (see for example, Brennan & Schwartz (1979: 152)) is: 3.5.1.2 The following notation is used in describing the FD methods: equal steps for some Δ t . The FD method estimate of ( ) 3.5.2 implEmEntation of tHE fD mEtHods 3.5.2.1 The following five paragraphs will summarise each of the five FD methods used in this paper. We refer the reader to eg to Wilmott (2006Wilmott ( : 1260 and references therein for a more elaborate treatment. 3.5.2.2 The explicit method involves discretising eq. (6) using the approximations in Appendix B4 as: The value of t ij U is the only unknown in eq. (7) and subsequently easily found, detailed in Appendix B5. Solving t ij U for all 0 < i, j ≤ n defines the new grid for each iteration. For the scheme to work, we must check that the 'probabilities' given by the parameters in eq. (6) connecting t ij U with the three known values ( ) 1 , The implicit method involves discretising eq. (6) using the approximations in Appendix B4 as: n 2 equations is obtained using eq. (8) for 0 < i, j ≤ n. The equations are solved simultaneously to find t ij U , 0 < i, j ≤ n in each iteration, with more detail in Appendix B6.
3.5.2.4 The discretised version of eq. (6) used in implementing the CN method is an average of eq. (7) and eq. (8). As with the implicit method, n 2 equations containing t ij U , with more detail in Appendix B7. 3.5.2.5 When using the Hopscotch method, each iteration consists of two stages. In the first stage, t ij U is calculated as was done in the explicit method using eq. (7) for all the odd or even points (a point is odd or even depending on whether its grid reference i j + is odd or even). The remaining grid points are then calculated in stage two as was done in the explicit method using eq. (7), but where the input grid is updated with the values from stage one (i.e. U +∆ − are replaced by the values calculated in stage one). The odd points are calculated in stage one and even points in stage two in the first iteration. The stages in which the odd and even points are calculated swap for each subsequent iteration. More detail is found in Appendix B8.
3.5.2.6 The ADI method also consists of two stages in each iteration. The first stage calculates a new grid at 1 2 t t + ∆ by solving θ 1 implicitly and θ 2 explicitly. The values of for 0 < i, j < n are found using the grid containing t t ij U +∆ for 0 < i, j ≤ n. The second stage calculates a new grid at time t by solving θ 2 implicitly and θ 1 explicitly, using the grid from stage one as an input. Each stage results in n independent sets of n simultaneous equations, the solution of which is the new grid. More detail is found in Appendix B9.

Comparisons of tHE fd mEtHods
3.5.3.1 Finite difference methods should preferably be stable, consistent and convergent. An FD method is stable if an error arising in any iteration remains bounded. An FD method is consistent if the truncation error ( ) 3.6 ODE methods 3.6.1 Five ODE methods used to estimate the solution to eq. (3) using eq. (4) are considered. Comparison of additional methods such as those proposed by Bulirsch & Stoer (1966) and Bashforth & Adams (1883) is an avenue for further research. The ODE methods involve iteratively calculating the values of ( ) ( ) and ending when T t τ = − (Shampine, 1994). The bond price is estimated as ( ) 3.6.2 An explicit method uses ( ) ( ) considered is the explicit Euler (or Euler/explicit) method. The implicit method considered is an implicit implementation of the Euler method and the CN method is an average of the explicit and implicit Euler methods. The Taylor method is explicit, although implicit and CN versions can be developed. The Runge-Kutta method considered is explicit, with implicit examples found in Alexander (1977).

3.6.3
The following functions of C 1 and C 2 are used for the ODE methods: 3.6.4 For the Euler method we use a first order Taylor expansion in each iteration as follows: 3.6.5 For the second order Taylor method we use a second order expansion in each iteration as follows: 3.6.6 For the fourth order explicit Runge-Kutta method each iteration is calculated as: 3.6.7 The implicit method requires ( ) τ + ∆ to be solved for simultaneously using the two equations: The value for ( ) A τ τ + ∆ in each iteration is then calculated as: 3.6.8 The CN method requires τ + ∆ to be solved for simultaneously using the two equations: The value for ( ) A τ τ + ∆ in each iteration is then calculated as: 3.6.9 The implicit and CN methods are expected to be more stable for larger τ ∆ than the explicit method (Granville, 1988). The local truncation error, defined as measures accuracy for each iteration. As 0 t ∆ → , LTE converges to 0 slower for the Taylor series than the Runge-Kutta method and even slower for the Euler method (Shampine, 1994).

METHODOLOGY
Without loss of generality, the numerical methods are compared for the sets of parameters in Table 1 below and the six times to maturity (in years) τ = 2, 5, 10, 15, 20, 30. To stress test the results, further comparisons are made for selected numerical methods across a broader random selection of fifty sets of parameters at all times to maturity τ = 1, 2, …, 30.

4.1
Selection of parameters 4.1.1 Five sets of parameters were selected to ensure generality across different yield curve shapes and different average bond yields. Increasing, decreasing and humped yield curves are represented within these sets of parameters, as well as average bond yields ranging from 0,02 i = p.a. to 0,13 i = p.a. The yield curves obtained from the five parameter sets are found in Appendix A.
4.1.2 Three of the five sets of parameters were chosen to roughly represent yield curves based on data obtained from Bloomberg for government bonds in the countries of South Africa, Brazil and the United States of America (USA) on 1 June 2016. Since the focus of this paper is not parameter estimation, the parameters are chosen only to roughly represent the market yield curves. There exists a substantial amount of literature on estimation procedures that could be used for more accurate estimates, including the papers by Chen & Scott (2003) and Pearson & Sun (1994).

4.1.3
We posit that the choice of the three countries was motivated by the fact that South Africa and Brazil are developing BRICS countries while the USA is a developed country. What is evident in all cases, as expected, is that zero-coupon bond prices decrease with term. 4.1.4 The other two sets of parameters were chosen to obtain a humped yield curve (the arbitrary set of parameters) and a decreasing yield curve (the check set of parameters). 4.1.5 The majority of the analysis is based on the five sets of parameters specified in Table 1. As previously mentioned, further analysis is done for selected numerical methods at all integer times to maturity τ = 1, 2, …, 30 based on an additional fifty sets of parameters. The additional sets of parameters are chosen randomly and used only for the methods appearing computationally efficient for the parameters in Table 1. The use of random parameters is based on the method used by Broadie & Detemple (1996), with these additional random parameters being selected as follows: The values for δ 0 , δ 1 and δ 2 are selected from a uniform distribution with range ( ) 0;0,04 randomly and independently . The values for μ 1 , μ 2 , λ 11 and λ 22 are selected randomly and independently from a uniform distribution with range ( ) 0;0,5 . The values for λ 12 and λ 21 are selected randomly and independently from a uniform distribution with range ( ) 0,5;0 − . Y 1 and Y 2 are set equal to 1. The random parameter selection is used to stress test the results and ensure generality of the results. The numerical methods are compared in this paper based on the two criteria of accuracy and computational efficiency.

4.2.2
The price of bonds considered in this paper assume a face value of USD1. Accuracy, however, is measured in dollars and cents for a bond with a face value of $100 (i.e. accuracy is measured to four decimal places for a bond with face value of USD1). Accuracy error is defined as

P t T is the true bond price and ( )
, P t T is the estimated price. For example, the accuracy for a bond with price of $0,8575 is measured assuming the face value is $100 and the price is $85,75. Two additional measures of accuracy are further used in this paper, with these being the root mean square relative error (RMS) and the mean absolute error (MAE). The RMS is defined as − . The RMS penalises large errors more than the MAE, and defines the error relative to the true bond price (Broadie & Detemple, 1996).

4.2.3
Computational efficiency is measured as the number of seconds taken to obtain the estimated bond price. Computational efficiency depends largely on the efficiency of the computer code used in implementing the numerical method, and is also influenced by the computer used to run the code. An effort was made to ensure the efficiency of the code, but the results may be impacted by inefficiencies in the code.

4.3
Determining a baseline for comparison of accuracy 4.3.1 The absence of a closed-form solution results in complexity in determining a true value to be used as a baseline for comparison of accuracy. The restriction of 12 21 0 λ λ = = in eq. (4), however, results in a system of ODEs for which a closed-form solution can be found. The check set of parameters was chosen with this restriction (last column of parameters), for which the closed-form solution (as calculated using WolframAlpha (2016)) is ( ) The five ODE methods all converge to the true solution (based on accuracy to seven decimal places) for the check set of parameters. For the remaining four sets of parameters and all additional random sets of parameters for which λ 12 <0 and λ 21 <0, the value to which the ODE methods converge (using sufficiently small Δ τ ) is used as the true value. It is further shown in the results that the other numerical methods converge to this value.

RESULTS
The results begin with comparisons across the three categories of methods, followed by detailed comparisons of the methods within each category.

5.1
Overview 5.1.1 The accuracy and computational efficiency of the numerical methods used to estimate the bond prices in Table 2 is summarised in Table 3. RMS and MAE are used as accuracy measures, calculated across the estimated bond prices for the six times to maturity for each parameter set and method. The total computational time required to estimate the six prices for each parameter set and method is used as the measure of efficiency.
5.1.2 The ODE methods have higher levels of accuracy while maintaining quicker computational times than the FD and simulation methods. The least efficient ODE method is a minimum of 14,7 times more computationally efficient than the most efficient of the FD and simulation methods over the five sets of parameters, while maintaining an RMS of 0 (rounded to four decimal places -referred to as 0 in the remainder of this paper). The simulation and FD methods are unable to achieve this level of accuracy in a reasonable computation time (explored further below). The ODE methods are clearly preferred on the basis of accuracy and efficiency but are only applicable to affine short-rate models.  iii. 5000 simulations for antithetic variates due to its lower variance; iv. ∆ t = 0,005 for simulation methods to achieve some degree of accuracy; v. ∆ t for the FD methods is chosen such that using ½∆ t decreased average e i by less than 0,01. Specifically, ∆ t equals: 0,01 for Hopscotch and implicit, 0,005 for explicit, 0,1 for ADI and 0,25 for CN; and vi. ∆ τ for the ODE methods is chosen as the maximum for which RMS = 0. The FD methods have higher levels of accuracy and quicker computational times than the simulation methods. The most efficient FD method is a minimum of 49,7 times more computationally efficient than the most efficient simulation method, while having a 94,8 per cent lower MAE across the five parameter sets and 50,8 per cent lower MAE on average. The FD methods have greater consistency in accuracy across the five parameter sets. These results are consistent with the comparison between FD methods and simulation methods found in Boyle (1977) who proposes that simulation methods are computationally inefficient. These results are also consistent with those in Wilmott (2006) in which it is argued that for models containing less than four random factors, finite difference schemes are more efficient than simulation methods. 5.1.4 The Euler, CN, Taylor and Runge-Kutta ODE methods as well as the ADI and CN finite difference schemes (selected due to their apparent computational efficiency) are analysed further using fifty additional randomly selected sets of parameters as previously described. The accuracy and computational efficiency of these methods when applied to the randomly selected sets of parameters is summarised and compared in Table 4, with the full results available in Appendix C. 5.1.4 The relative ranking of theses numerical methods based on efficiency and accuracy does not differ between the parameter sets in Table 1 and the random parameter sets. The conclusions drawn above based on the results for the sets of parameters in Table 1 hold when considering the random selected sets of parameters. In particular, the ODE methods remain more efficient and accurate across the random parameter sets. In Table 4 it can be seen that the least efficient of the four ODE methods shown is more efficient than the most efficient FD method on average, while the most efficient ODE method is more than 400 times more efficient on average. Furthermore, the least accurate ODE method has a 98 per cent lower MAE than the most accurate FD method on average. An improvement in the accuracy of the FD methods causes a significant worsening of their already low relative efficiency. Runge 0,0000 0,0000 0,00 0,0000 0,0000 0,00 0,0000 0,0000 0,00 0,0000 0,0000 0,00 0,0000 0,0000 0,00 Euler 0,0000 0,0000 0,14 0,0000 0,0000 0,06 0,0000 0,0000 0,06 0,0000 0,0000 0,08 0,0000 0,0000 0,05 Taylor 0,0000 0,0000 0,22 0,0000 0,0000 2,56 0,0000 0,0000 0,05 0,0000 0,0000 0,47 0,0000 0,0000 0,01 Implicit ( Tables 3 and 4, the Runge-Kutta method is the most computationally efficient method for all sets of parameters considered. Table 5 further considers the relative efficiency of the ODE methods when requiring a greater degree of accuracy. 2 i. Grid size = 1600 for the two FD methods; ii. ∆ t for the two FD methods is the largest value for which a reduction in ∆ t by a factor of two results in a reduction in e i by less than 0,01; and iii. ∆ τ for the ODE methods is the maximum value for which e i < 0,001.

5.2.2
The Runge-Kutta method estimates bond prices accurately (measured to six decimal places) in a computation time of less than 0,01 seconds for all 1 30 τ ≤ ≤ for all parameter sets in Table 1, on average 10 times quicker than the second most efficient method (determined separately for each parameter set). The Runge-Kutta method further estimates bond prices accurately (e i < 0,001) for all times 1 30 τ ≤ ≤ for 50 random parameter sets within 1.5 seconds. The other four ODE methods are compared in Figure 1 in addition to the results above by considering the rates of convergence to the true bond price (Runga-Kutta is excluded due to its efficiency making graphical comparison less relevant).

5.2.3
For each of the five parameter sets in Table 1, the relative efficiencies of the methods are consistent whether computation time is based on the time required to achieve an RMS of 0 or MAE of less than 0,0000005 (Table 3 and Table 5 respectively). The CN method converges quicker than the Euler method which in turn converges quicker than the implicit method (for all five parameter sets). The Euler method requires between two and six per cent of the computation time required by the implicit method. The efficiency of the CN method relative to the Euler method is less consistent, with the CN method requiring between two and 75 per cent of the computation time required by the Euler method. The efficiency of the Taylor method depends on the particular set of parameters. The Taylor method appears to be more efficient (relative to the efficiency of the other methods) if λ 12 and λ 21 are close to 0, This is the case for the USA and check parameter sets, where the Taylor method is more efficient than the Euler method. On average, the CN method converges the second quickest of the ODE methods, although the Taylor method converges marginally quicker for the check parameter set. The efficiency of the CN method is due to the target level of accuracy being obtained using a larger Δ τ than the explicit and Taylor methods, which appears to offset the complexity of needing to solve two simultaneous equations. 5.2.4 Similar results are obtained for the random parameter sets (the implicit method is excluded). The Runge-Kutta method is on average 8,9 times more efficient than the second most efficient method (determined independently for each parameter set). The CN method is the second most efficient method for 84 per cent of the parameter sets, and is more efficient than the Euler method for 88 per cent of parameter sets. The CN method, on average, is 5,1 times more efficient than the Euler method and 17,8 times more efficient than the Taylor method. The relative ranking of the Taylor and Euler methods is dependent on the parameter set. On average, however, the Euler method is more efficient than the Taylor method, being 6,9 times more efficient on average and more efficient for 72 per cent of the parameter sets.

5.3
Finite difference scheme methods 5.3.1 Finite difference scheme methods can be adapted to non-affine models of the short rate, unlike the ODE methods which rely on the affine structure of the model considered in this paper (Shreve, 2004). We do not suggest the results in this paper will hold for non-affine models without having tested this assertion. We merely propose that the results may be of interest in spite of having previously shown the superior efficiency of the ODE methods. We further consider the results in comparison to the previous studies mentioned in ¶3.3. 5.3.2 The accuracy and efficiency of the five FD methods depends on both the grid size and Δ t . An understanding of the accuracy of the FD methods is first considered in Table 6, after which the efficiency to achieve similar levels of accuracy is considered. Table 6 represents how the accuracy of these methods changes for selected sets of parameters with an increase in the size of the grid and reduction in the size of the time step used Δ t . 5.3.3 As can be seen in Table 6, assuming the same grid size is used that is not too small, the five FD methods converge to the same estimated bond price as Δ t is reduced (referred to as the limited estimate in this paper), but the rate of convergence differs between the methods. The implicit and explicit methods converge to the limited estimate from opposite directions. The size of the grid determines the level of accuracy attainable with a reduction in Δ t (i.e. the accuracy of the limited estimate). An increase in the size of the grid increases the accuracy of the limited estimate. 5.3.4 Given the five FD methods converge to the same estimated bond price for the same size of grid, the efficiency of the FD methods can be compared for a defined grid size. This is done in the following paragraphs. Although the relationship between grid size and accuracy is an important consideration in overall efficiency of the FD methods, it is not a core result of this paper apart from what is discussed in ¶5.3.10 and as such is explored in more detail in Appendix D.
5.3.5 Table 7 compares the efficiency of the FD methods for a defined grid size by comparing the computational time required to achieve a given level of accuracy through a reduction in Δ t . The Δ t required to achieve the defined level of accuracy is also shown. 5.3.6 The results in Table 7 are analysed in the paragraphs below, but the distinction between implicit 1 and implicit 2 is first considered as this impacts the results. For all parameter sets and times to maturity shown in Table 7, the implicit method has an estimated bond price that decreases towards the true bond price with a reduction in Δ t but then away from the true bond price towards the limited estimate. The results labelled implicit 1 are for Δ t being selected as the maximum value for which e i is less than 0,0025 greater than that for the limited estimate. The results labelled implicit 2 are for Δ t being selected as the maximum value for which the difference between the estimated bond price and the limited estimate is less than 0,0001. Implicit 2 is preferred for comparison as the limited estimate is the estimated price to which the FD methods converge, and as such considered the base value around which comparisons of efficiency of the FD methods should be made. Comparisons of the implicit method in the following paragraphs refer to implicit 2 .  Table 7 shows that the ADI method is the most computationally efficient of the FD methods across all five parameter sets, followed by the CN method. The Hopscotch method is more efficient than the explicit method. The implicit method is the least efficient method on average, with its theoretical stability not seeming to translate into improved computational efficiency. The time step required to obtain a given level of accuracy (relative to the limited estimate) is the largest for the CN and ADI methods, significantly larger than that required by the other three methods (based on the Δ t values used for Table 7). This offsets the complexity of these two methods (although the ADI method is significantly less complex), resulting in their efficiency. 5.3.8 Based on the computation times in Table 7, the explicit method requires on average 38 per cent of the computation time required by the implicit method, while the Hopscotch method requires on average 33 per cent of the computation time required by the explicit method. The ADI method requires on average 3.5 per cent of the computation time required by the Hopscotch method and 15 per cent of that required by the CN method. For the random parameter sets (Table 4), the ADI method is six times more efficient on average than the CN method with similar levels of accuracy. 5.3.9 As detailed in ¶3.3, Geske & Shastri (1985) showed the explicit method to take roughly 60 per cent of the computation time taken by the implicit method, while Hull & White (1990) argue that the explicit scheme uses between 40 and 70 per cent of the computation time of the implicit method. The results in this paper similarly show the explicit method to be computationally more efficient. The results in Table 7 show the explicit method requiring between 25 and 75 per cent of the computation time of the implicit method. Previous results when using these methods (see eg Cairns (2004)) further showed that the CN method converges more quickly than the explicit and implicit method, which is the result of this paper.
5.3.10 The analysis in ¶ ¶5.3.7 and 5.3.8 has shown the ADI method to be the most efficient of the FD methods. The efficiency of this method is considered in more detail in Table 8 to draw a comparison relative to the efficiency of the Runga-Kutta method (the most efficient ODE method). The efficiency is considered for the Brazil parameter set as the smallest grid size was required to achieve a given degree of accuracy for this set of parameters.  The ADI method is unable to obtain an RMS of 0 within one hour of computation time. This is in contrast to the Runge-Kutta method which obtains an RMS of 0 within 0.01 seconds of computation time. Given the ADI is the most efficient of the FD methods, it can be concluded that the FD methods are unable to obtain the accuracy of the ODE methods (in this case measured as an RMS of 0) within a reasonable computation time.

5.4
Simulation methods 5.4.1 This section considers the simulation methods in more detail due to their flexibility. In particular consideration is given to the impact of the number of simulations and the size of the Δ t on the accuracy and efficiency of the simulation methods. As mentioned previously, these impact the two causes of error in the estimated price ( ) , P t T for the simulation methods. The number of simulations will determine the variation in the sample estimate of the expectation. Δ t will determine the accuracy of the estimate as the discretisation error of the short rate and the error in estimating an integral with a discrete sum reduces as Δ t → 0. The variation in the sample estimate as measured with standard deviation is analysed first. This is done by considering the change in the standard deviation of the estimates with a change in the number of simulations for a fixed Δ t , enabling a comparison of the relative efficiency of the two simulation methods. The impact of Δ t on the accuracy of the simulation methods is then considered in ¶ ¶ 5.4.5 and 5.4.6. 5.4.2 Table 9 below compares the reduction in the observed and theoretical standard deviations of the estimated bond price with an increase in the number of simulations the South African parameter set.  Table 9 were found when considering the Brazil (setting Δ t = 10) and Check (setting Δ t = 5) parameter sets, with the standard deviation of the estimated bond price being between 40 and 70 per cent lower using the antithetic variates method as opposed to Monte Carlo simulation, assuming the same Δ t and number of simulations. Furthermore, the observed sample standard deviation does not differ significantly from the theoretical standard deviation of the estimated bond price in any of the parameter sets considered. The theoretical standard deviation is a quick measure of the number of simulations required to reduce the sample variance to within an acceptable limit, and highlights the inefficiency if very low standard deviation is required. that required by Monte Carlo simulation by a factor of 2. It, therefore, has an estimated 28 per cent lower computation time (based on assuming that it results in a 40 per cent lower variance than Monte Carlo simulation). Similarly, if the antithetic variates method results in a 70 per cent lower variance than Monte Carlo simulation, it will have an estimated 82 per cent lower computation time. It, therefore, has a computation time of between 28 and 80 per cent lower than that of Monte Carlo simulation in order to obtain the same level of variance in the estimated bond price. This is in spite of the antithetic variates method requiring approximately double the computation time of Monte Carlo simulation for the same Δt and number of simulations. This is as a result of the lower variance offsetting the doubling of the computation time. 5.4.5 The above has analysed the effect of increasing the number of simulations on the standard deviation of the estimate. For a given number of simulations, the accuracy of the estimates will depend on Δ t . This is considered in Table 10. The error decreases as the time step decreases. The level of accuracy, however, depends on the set of parameters considered. This is similar to the results in Table 2, where the average accuracy (using MAE) across the South African, Brazil and USA parameter sets is more than 10 times higher than that across the check and arbitrary parameter sets using the same Δ t and number of simulations.

CONCLUSIONS
6.1 Numerical methods used to estimate the price of zero-coupon bonds in a two-factor CIR model were compared on the basis of accuracy and computational efficiency. Five ODE methods, five FD methods and two simulation methods were considered. The ODE methods were shown to be significantly more efficient, consistent and accurate than the FD methods, which in turn were more efficient, consistent and accurate than the simulation methods. The ODE methods can only be used in affine short-rate models, while the FD methods can be applied to many non-affine short-rate models (Shreve, 2004).

6.2
The Runge-Kutta method is the most preferred of the numerical methods considered, having the quickest computation time and highest accuracy for all parameters and times to maturity considered. The Runge-Kutta method is, on average, between nine and ten times more computationally efficient than the second most efficient method. This paper would suggest that the lower theoretical truncation error of the Runge-Kutta method translates into improved efficiency. The Crank-Nicolson ODE method is the second most computationally efficient method, while the relative efficiency of the Taylor and Euler methods is dependent on the parameter set considered. The implicit method is the least efficient of the ODE methods, in spite of its theoretical stability.

6.3
The FD methods are unable to obtain the same level of accuracy as the ODE methods in a reasonable computation time. This is the result of the accuracy of the FD methods being largely constrained by the size of grid used. The ADI method is the most efficient FD method and was on average six times more efficient than the Crank-Nicolson method, the second most efficient FD method. The use of the antithetic variates methods improved the efficiency of Monte Carlo simulation, but these two methods were the most inefficient.

6.4
The numerical methods used in this paper are some of the most commonly used numerical methods in mathematics of finance. In mathematics of finance derivative pricing, a bigger proportion of derivatives are over the counter traded and their prices sometimes may not be obtained in closed-form due to their exotic nature. Our paper contributes to the literature in terms of providing evidence that the chosen numerical methods, though applied to a particular bond pricing problem, could be used in cases where closed-form solutions either do not exist or are difficult to obtain using classical methods. We demonstrated the advantages of each method based on efficiency and accuracy. We acknowledge that recent numerical methods, eg Fast Fourier transform (FFT), Cosine method (COS) etc, require knowledge of the distribution or characteristic function of the underlying, which is one of their weaknesses.

APPENDIX A Yield curves
The following diagrams illustrate the shape of the yield curves arising from the five sets of parameters found in Table 1.

B1.
Description of how eq. (6) is obtained B1.1 The transformation in eq. (5), as presented by Brennan & Schwartz (1979: 152), results in the following: The partial derivatives of f with respect to Y 1 and Y 2 as functions of θ 1 and θ 2 are obtained as follows: 3 It is now possible to specify eq. (2) in terms of θ 1 and θ 2 as follows: Eq. (6) is subsequently obtained using a 0 , a 1 , a 2 , b 1 and b 2 for notational compactness.

B4.
Approximations used in discretisation of eq. (6) for the FD methods The following approximations are used in discretising eq. (6) to obtain eq. (7) and eq. (8). Backward differences are used for

B5.
Solving for t ij U in the explicit method B5.1 The explicit method is developed from Cairns (2004:169). The value of t ij U for 0 < i, j < n is found using eq. (7) as follows: .
B7.2 Eq. (11) for 0 , i j n < < can be represented as: . Q 1 and Q 2 are n × n co-efficient matrices obtained from eq. (12) and eq. (13) above. X 1 and X 1 are n × 1 target matrices obtained from eq. (12) and eq. (13) above. Each row and column index x in Q 1 and each row index x in X 1 specifies the value relating to the ( ) , i j grid point for fixed j and 0 i n < ≤ . Each row and column index x in Q 2 and each row index x in X 2 specifies the value relating to the ( ) , i j grid point for fixed i and 0 j n < ≤ . Accuracy and efficiency for the fifty random parameter sets for selected numerical methods. Accuracy is measured using RMS and MAE and efficiency is measured using total calculation time in seconds, calculated over thirty integer times to maturity 1,...,30 τ =

D1.
The overall efficiency of the FD methods is dependant on the grid size. An increase in the size of the grid increases the accuracy of the limited estimate but also increases the computational time of the method. The tables below consider the change in the computation time with a change in the size of the grid for selected sets of parameters at selected times to maturity. The multiplier represents the factor at which the computational time for a given size of grid increased relative to the previous grid size, i.e.  The computation time for the explicit, Hopscotch and ADI methods increase largely in direct proportion to an increase in the grid size as shown above in Tables 12 and 13. The computation time for the implicit and CN methods increases exponentially with an increase in grid size as a result of needing to create an n × n matrix to solve the n simultaneous equations, where n is the size of the grid. If an increase in the level of accuracy of the limited estimate is required (and subsequently a larger grid size needs to be used), the CN and implicit methods