Homotopy perturbation transform method for pricing under pure diffusion models with affine coefficients

1 Introduction

Stochastic processes have taken over the world of financial modelling. Starting with simple Geometric Brownian Motion well described by Bachelier (1900), to more sophisticated processes for better fitness and calibration of market fluctuations. A huge range of papers have considered Lévy processes as their driving force, they usually take the umbrella of jump diffusion process. Basically the state process $X_t$ follows a Brownian motion with drift for which a pure jump process, usually of Poisson type, is added on in order to accommodate abrupt changes in the market. In this work we restrict ourselves to pure diffusion processes for which the drift, the variance as well as the interest rate processes are all affine functions of $X_t$ ; we refer the reader to Eq. (2.4) in the next section for the explicit form. Duffie et al. (2000) present a general framework of diffusion processes under affine coefficients (also termed as affine models) then apply it to a two dimensional option pricing problem. These diffusion models have the advantage of providing tractability of closed form formula for a wide range of asset price such as fixed income securities: bonds, options and swaps. Technically, in dealing with these diffusion models we apply transforms that will result later in a system of ordinary differential equations of Riccati type that can be solved analytically or numerically in order to compute the asset price. In most cases, numerical methods using Fourier transform or inverse Laplace transform are needed to solve the integro-differential equation that arise thereof.

Various types of diffusion models have been introduced in the literature and used in different sectors. They differ in the model parameter while preserving some properties that are essential in asset price dynamics such as mean reversion.

In one-factor models the volatility process is a deterministic function of t, this includes the class of affine term structures that is encountered in interest rates and option pricing contexts, see Duffie (2005), Björk (2004), Black and Scholes (1973) and Merton (1974).

Heston (1993) introduces a two-factor model where the volatility is now stochastic, to accommodate the implied volatility smile encountered in the financial markets as shown by Kotzé et al. (2015). Bates (1996) extends the Heston model in adding a up jumps into the diffusion process to model a huge flux of information that can occur in the market.

Fang (2000) brings a three-factor model with jump diffusion process that looks quite stable. The model takes its roots from the Bates model to which an extra parameter is added: the long-term volatility, which is important for instruments like bonds that have long term maturity. Banks are often interested in this parameter. Christoffersen et al. (2009) consider a three-factor model with no jumps, but the state process is driven by a deterministic drift, plus two Brownian motions. However, as the number of factors increases one must expect the model to become more robust, but less realistic. In this paper we present a detailed framework of a class of asset prices whose pay-off at a future time (the maturity time) T is of the form $e^{u{X}_T}$ where $X_t$ is of a pure diffusion type with affine coefficients. This class includes fixed income security prices. A quick overview of the pure diffusion with affine coefficients is provided together with an application on bond price in the context of single, two and three-factor models based on the Heston type. In other words, the volatility is considered to be stochastic. This Gives them a wider range of application. Bravo (2008) uses an affine process for pricing longevity bonds; Jang (2007) applies a two-factor affine process in insurance; Crosby (2008) uses it in pricing a class of exotic commodities and options in a multi-factor model.

In addition we introduce a modified form of homotopy perturbation and Laplace transform methods to value financial models of diffusion type with affine coefficients. These models lead to the need of solving systems of Riccati differential equations. Laplace transform method (LTM) alone is incapable of handling such equations, instead some variants of the LTM prove to better handle nonlinear differential equations. Among those variants we cite the Laplace decomposition algorithm (Khuri, 2001; Khan, 2009) and the H2LTM, (Fatoorehchi and Abolghasemi, 2016) which are obtained with the help of the Adomian decomposition and Adomian polynomials respectively. Another successful variant of LTM is obtained in coupling it with variation iteration methods. Alawad et al. (2013) used it to solve space–time fractional telegraphic equation as it allowed them to overcome the difficulty arising from finding Lagrange multipliers. LTM has also known great success when combined with differential transform methods on solving non–homogenous equations, see Alquran et al. (2012). On the other hand, the homotopy perturbation method is a combination of the classical perturbation technique and the homotopy technique whose origin is in topology , more on homotopy can be found in Hilton (1953), but not restricted to small parameters as it occurs with traditional perturbation methods. For example, the HPM requires neither small parameters nor linearisation, but only few iterations to obtain highly accurate solutions. The standard homotopy perturbation method was proposed by He (1999) as a powerful tool to approach various kinds of nonlinear problems. It can also be viewed as a special case of Homotopy Analysis Method (HAM) proposed by Liao (1992, 1997). For the past decade, many improvements on HAM have been introduced, one of it being the Homotopy Analysis Transform Method (HATM) which is basically a HAM coupled with a Laplace transform. The method is very powerful, fast converging and accurate. Recently, it has been applied in many different areas of science including fluid dynamics, wave theory (Kumar et al., 2014a, 2015b, 2016b), quantum physics, see Kumar (2014), and many more, with extension to fractional cases. Kumar et al. (2014c) applied this method on Volterra integral equation to obtain good quality results. Another improvement of HAM is obtained by coupling it with the Samudu transform which gives rise to the Homotopy Analysis Samudu Transform Method (HASTM), see Kumar and Sharma (2016); Kumar et al., 2016a for more details. Likewise, the Samudu transform has also been introduced in HPM to generate the Homotopy Perturbation Samudu Transform Method (HPSTM), see Singh et al. (2013); Singh et al., 2014b. Singh et al. (2014a) used the method successfully to get analytical and numerical solutions of nonlinear fractional equations found in the area of biological population model. Also, Kamdem (2014) proposed a generalised integral transform based on the homotopy perturbation method where various integral transforms were used. In this paper we are interested in the combination of the HPM with Laplace transform giving rise to the Homotopy Perturbation Transform Method (HPTM). The method has shown success already in obtaining solutions of the Navier–Stokes equations (Kumar et al., 2015a), gas dynamics equations coming from fluid dynamics in the case of fractional differential equation as explored by Kumar et al. (2012). The method has also shown success in solving KdV equations arising in wave theory (Goswami et al., 2016) as well as Fokker-Planck equations commonly found in solid-state physics (Kumar, 2013). Another useful application of the technique is found in Kumar et al. (2014b) in which the authors derive the price of a plain vanilla call option of European type under Black–Scholes model in the financial market.

This paper is structured as follows: Section 2 reviews the formalism of diffusion models with emphasis on the case with affine coefficients. In Section 3, we introduce the basic concept of homotopy perturbation transform method (HPTM). In Section 4, we describe the solution procedure of the HPTM for interest rate models, especially the two and three-factor stochastic volatility models. Finally, the conclusions are presented in Section 5.

2 Mathematical description of affine models

Consider the financial market model $\mathcal{M}=\left(\mathrm{\Omega }\text{,}\mathcal{F}\text{,}\mathbb{P}\text{,}{\left({\mathcal{F}}_t\right)}_{t⩾0}\text{,}{\left(S_t\right)}_{t⩾0}\right)$ where $\mathrm{\Omega }$ is the set of all possible outcomes of the experiment known as the sample space, $\mathcal{F}$ is the set of all events, i.e. permissible combinations of outcomes, $\mathbb{P}$ is a map $\mathcal{F}⟶[0\text{,}1]$ which assigns a probability to each event, ${({\mathcal{F}}_t)}_{t⩾0}$ is a natural filtration and $S_t$ a risky underlying asset price process. The triplet $(\mathrm{\Omega }\text{,}\mathcal{F}\text{,}\mathbb{P})$ is defined as a probability space. Let ${(W_t)}_{t⩾0}$ denote a $\mathbb{P}$ -Wiener process, $\sigma >0$ the volatility of the underlying asset, $\mu (t\text{,}{X}_t)$ the drift parameter. Suppose $X_t$ satisfies the following stochastic differential equation

Under an equivalent martingale measure $\mathbb{Q}$ , the price $\psi (t\text{,}x)$ at time t of a contingent claim that pays off ${\mathrm{\Phi }}_T$ at maturity time $T⩾t$ is given by

Let us consider an auxiliary process $\mathrm{\Psi }(t\text{,}{X}_t)=E^Q\left(e^{-{\int }_0^{T}r(s\text{,}{X}_s)\mathit{ds}}{\mathrm{\Phi }}_T|{\mathcal{F}}_t\right)$ . For simplicity of notation, we may from time to time write ${\mathrm{\Psi }}_t$ to denote $\mathrm{\Psi }(t\text{,}{X}_t)$ . Same for other variables as well. Applying the Ito’s differentiation we get $$\psi (t\text{,}x)-\psi (0\text{,}x)={\int }_0^t\frac{\partial {\mathrm{\Psi }}_s}{\partial s}\mathit{ds}+{\int }_0^t\frac{\partial {\mathrm{\Psi }}_s}{\partial X_s}{\mathit{dX}}_s+\frac{1}{2}{\int }_0^t\frac{{\partial }^2{\mathrm{\Psi }}_s}{\partial X_s^2}d{[X\text{,}X]}^c={\int }_0^t\frac{\partial {\mathrm{\Psi }}_s}{\partial s}\mathit{ds}+{\mu }_s\frac{\partial {\mathrm{\Psi }}_s}{\partial X_s}+\frac{1}{2}{\sigma }_s{\sigma }_s^{\top }\frac{{\partial }^2{\mathrm{\Psi }}_s}{\partial X_s^2}\mathit{ds}+{\int }_0^t{\sigma }_s\frac{\partial {\mathrm{\Psi }}_s}{\partial X_s}{\mathit{dW}}_s={\int }_0^t\frac{\partial {\mathrm{\Psi }}_s}{\partial s}\mathit{ds}+{\mu }_s\frac{\partial {\mathrm{\Psi }}_s}{\partial X_s}+\frac{1}{2}{\sigma }_s{\sigma }_s^{\top }\frac{{\partial }^2{\mathrm{\Psi }}_s}{\partial X_s^2}\mathit{ds}+{\int }_0^t{\sigma }_s\frac{\partial {\mathrm{\Psi }}_s}{\partial X_s}{\mathit{dW}}_s={\int }_0^t\frac{\partial {\mathrm{\Psi }}_s}{\partial s}\mathit{ds}+{\mu }_s\frac{\partial {\mathrm{\Psi }}_s}{\partial X_s}+\frac{1}{2}{\sigma }_s{\sigma }_s^{\top }(s\text{,}{X}_t)\frac{{\partial }^2{\mathrm{\Psi }}_s}{\partial X_s^2}\mathit{ds}+{\int }_0^t{\sigma }_s\frac{\partial {\mathrm{\Psi }}_s}{\partial X_s}{\mathit{dW}}_s$$

Under no-arbitrage conditions, the discounted pay-off ${\mathrm{\Psi }}_t$ , must be martingale, meaning the drift part must be zero. That is,

In Affine framework as described by Cont and Tankov (2004), we consider $\mu \text{,}\sigma \text{,}\mathrm{and}r$ to be Affine in X. That is,

Eq. (2.7) together with initial condition (2.8) is a nonlinear system of ODE of Riccati type. In general, Riccati equations do not have analytical solutions, hence numerical methods have to be used. In addition, these equations have been reported to be stiff. Hence the use of explicit methods will require a high mesh refinement to produce acceptable solutions. This will result in an increase in the computational cost. In this article we propose a Laplace Transform Homotopy Perturbation method to circumvent the stiffness problem.

3 Basic idea of homotopy perturbation transform method

We first define the Laplace transform (LT) and its inverse transform, and list useful properties employed in this paper.

4 Numerical experiments

4.1

4.1 Interest rate models

We first look at affine term structure models that are found in interest rate models (see Björk, 2004). Here the state process is the interest rate itself r following the dynamics: $$\mathit{dr}(t)=\mu (t\text{,}r(t))\mathit{dt}+\sigma (t\text{,}r(t)){\mathit{dW}}_t$$ where $\mu (t\text{,}r(t))=a(t)r(t)+b(t)\text{,}{\sigma }^2(t\text{,}r(t))=c(t)r(t)+d(t)$ and $a(t)\text{,}b(t)\text{,}c(t)\text{,}d(t)$ are deterministic functions of t. This suggests that the state process $X_t$ corresponds exactly to the interest rate $r(t)$ . It is one dimensional, and as a result we have the following matching $$\begin{array}{c}\hfill K_0=b(t)\text{,}{H}_0=d(t)\text{,}{\rho }_0=0\text{,}\\ \hfill K_1=a(t)\text{,}{H}_1=c(t)\text{,}{\rho }_1=1\text{.}\end{array}$$

For a zero coupon bond that pays 1 at maturity T, we see that its price at any time $\tau$ prior to maturity is given by $$p(\tau \text{,}r)=e^{\alpha (\tau )+\beta (\tau )r(\tau )}\text{,}$$ with $\alpha$ and $\beta$ satisfying the system of stiff ODEs of the form

(4.1)

$$\left\{\begin{array}{c}\frac{\partial \alpha }{\partial \tau }(\tau )=b(\tau )\beta (\tau )+\frac{1}{2}d(\tau ){\beta }^2(\tau )\hfill \\ \frac{\partial \beta }{\partial \tau }(\tau )=-1+a(\tau )\beta (\tau )+\frac{1}{2}c(\tau ){\beta }^2(\tau )\hfill \end{array}\right.\text{,}$$ where $\tau =T-t$ and the initial conditions are given by $$\alpha (0)=0\mathrm{and}\beta (0)=u\text{.}$$

We investigate numerical solutions of the system (4.1) by means of the HPTM. To this end we construct the following homotopy

(4.2)

$$\left\{\begin{array}{c}{A}^{\prime }(\tau )-{\alpha }_0(\tau )+p\left[{\alpha }_0(\tau )-b(\tau )B(\tau )-\frac{1}{2}d(\tau )B^2(\tau )\right]=0\hfill \\ B^{\prime }(\tau )-{\beta }_0(\tau )+p\left[{\beta }_0(\tau )+1-a(\tau )B(\tau )-\frac{1}{2}c(\tau )B^2(\tau )\right]=0\hfill \end{array}\right.\text{.}$$

Applying the Laplace transform on both sides of (4.2), we have

(4.3)

$$\left\{\begin{array}{c}\mathcal{L}\left\{A^{\prime }(\tau )-{\alpha }_0(\tau )+p\left[{\alpha }_0(\tau )-b(\tau )B(\tau )-\frac{1}{2}d(\tau )B^2(\tau )\right]\right\}=0\hfill \\ \mathcal{L}\left\{B^{\prime }(\tau )-{\beta }_0(\tau )+p\left[{\beta }_0(\tau )+1-a(\tau )B(\tau )-\frac{1}{2}c(\tau )B^2(\tau )\right]\right\}=0\hfill \end{array}\right.\text{.}$$

Using the differential property of the Laplace transform we have

(4.4)

$$\left\{\begin{array}{c}s\mathcal{L}\{A(\tau )\}-A(0)=\mathcal{L}\left\{{\alpha }_0(\tau )-p\left[{\alpha }_0(\tau )-b(\tau )B(\tau )-\frac{1}{2}d(\tau )B^2(\tau )\right]\right\}\hfill \\ s\mathcal{L}\{B(\tau )\}-B(0)=\mathcal{L}\left\{{\beta }_0(t)-p\left[{\beta }_0(\tau )+1-a(\tau )B(\tau )-\frac{1}{2}c(\tau )B^2(\tau )\right]\right\}\hfill \end{array}\right.\text{.}$$

By applying inverse the Laplace transform on both sides of Eq. (4.4) and after algebraic simplification we get

(4.5)

$$\left\{\begin{array}{c}A(\tau )={\mathcal{L}}^{-1}\left\{\frac{1}{s}\left(A(0)+\mathcal{L}\left\{{\alpha }_0(\tau )-p\left[{\alpha }_0(\tau )-b(\tau )B(\tau )-\frac{1}{2}d(\tau )B^2(\tau )\right]\right\}\right)\right\}\hfill \\ B(\tau )={\mathcal{L}}^{-1}\left\{\frac{1}{s}\left(B(0)+\mathcal{L}\left\{{\beta }_0(t)-p\left[{\beta }_0(\tau )+1-a(\tau )B(\tau )-\frac{1}{2}c(\tau )B^2(\tau )\right]\right\}\right)\right\}\hfill \end{array}\right.\text{.}$$

Suppose the solution of Eq. (4.2) to have the following form

(4.6)

$$\left\{\begin{array}{c}A(\tau )=A_0(\tau )+{\mathit{pA}}_1(\tau )+p^2{A}_2(\tau )+\cdots \hfill \\ B(\tau )=B_0(t)+{\mathit{pB}}_1(\tau )+p^2{B}_2(\tau )+\cdots \hfill \end{array}\right.\text{,}$$ where $A_j(\tau )\text{,}{B}_j(\tau )\text{,}j=1\text{,}2\text{,}\dots$ are unknown functions which should be determined. Substituting Eq. (4.6) into Eq. (4.5), collecting the same powers of p and equating each coefficient of p to zero, results in

(4.7)

$$p^0:\left\{\begin{array}{c}{A}_0(\tau )={\mathcal{L}}^{-1}\left\{\frac{1}{s}\left(A(0)+\mathcal{L}\left\{{\alpha }_0(\tau )\right\}\right)\right\}\hfill \\ B_0(\tau )={\mathcal{L}}^{-1}\left\{\frac{1}{s}\left(B(0)+\mathcal{L}\left\{{\beta }_0(\tau )\right\}\right)\right\}\hfill \end{array}\right.$$

(4.8)

$$p^1:\left\{\begin{array}{c}{A}_1(\tau )={\mathcal{L}}^{-1}\left\{-\frac{1}{s}\mathcal{L}\left\{{\alpha }_0(\tau )-b(\tau )B_0(\tau )-\frac{1}{2}d(\tau )B_0^2(\tau )\right\}\right\}\hfill \\ B_1(\tau )={\mathcal{L}}^{-1}\left\{-\frac{1}{s}\mathcal{L}\left\{{\beta }_0(\tau )+1-a(\tau )B_0(\tau )-\frac{1}{2}c(\tau )B_0^2(\tau )\right\}\right\}\hfill \end{array}\right.$$

(4.9)

$$\begin{array}{cc}\hfill & ⋮\hfill \\ \hfill & p^j:\left\{\begin{array}{c}{A}_j(\tau )={\mathcal{L}}^{-1}\left\{-\frac{1}{s}\mathcal{L}\left\{-b(\tau )B_{j-1}(\tau )-\frac{1}{2}d(\tau ){\displaystyle \sum _{k=0}^{j-1}}{B}_k(\tau )B_{j-k-1}(\tau )\right\}\right\}\hfill \\ B_j(\tau )={\mathcal{L}}^{-1}\left\{-\frac{1}{s}\mathcal{L}\left\{-a(\tau )B_{j-1}(\tau )-\frac{1}{2}c(\tau ){\displaystyle \sum _{k=0}^{j-1}}{B}_k(\tau )B_{j-k-1}(\tau )\right\}\right\}\hfill \end{array}\right.\text{.}\hfill \end{array}$$

Example 4.1

We consider a particular case of the Vasicek model (see Björk, 2004). This model is obtained from Eq. (4.1) when the parameters are $b(t)=b\text{,}d(t)={\sigma }^2\text{,}a(t)=-a\text{,}c(t)=0$ and $u=0$ . The exact solution of the Vasicek model was found to be of the form

(4.10)

$$\beta (\tau )=\frac{1}{a}\left(e^{-a\tau }-1\right)\mathrm{and}\alpha (\tau )=\frac{(\beta (t)-\tau )\left(\mathit{ab}-\frac{1}{2}{\sigma }^2\right)}{a^2}-\frac{{\sigma }^2\beta (\tau )}{4a}\text{.}$$

The Taylor expansions of both $\alpha$ and $\beta$ at about zero at order 6 is $$\alpha (\tau )=-\frac{b}{2}{\tau }^2+\left(\frac{\mathit{ab}}{6}+\frac{d}{6}\right){\tau }^3-\left(\frac{a^{2}b}{24}+\frac{\mathit{ad}}{8}\right){\tau }^4+\left(\frac{a^{3}b}{120}+\frac{7a^2}{120}\right){\tau }^5-\left(\frac{a^{4}b}{720}+\frac{a^{3}d}{48}\right){\tau }^6+O{[\tau ]}^7\text{,}$$ and $$\beta (\tau )=-\tau +\frac{1}{2}a{\tau }^2-\frac{a^2}{6}{\tau }^3+\frac{a^3}{24}{\tau }^4-\frac{a^4}{120}{\tau }^5+\frac{a^5}{720}{\tau }^6+O{[\tau ]}^7\text{.}$$

To obtain the numerical solution of the (4.1) for the Vasicek model, we assume ${\alpha }_0(\tau )=A(0)=\alpha (0)=0$ and ${\beta }_0(\tau )=B(0)=\beta (0)=0$ . Solving Eqs. (4.7)–(4.9) for $A_j(\tau )\text{,}{B}_j(\tau )\text{,}j=0\text{,}1\text{,}\dots \text{,}6$ leads to the following result

(4.11)

$$A(\tau )=-\frac{b}{2}{\tau }^2+\frac{\mathit{ab}}{6}{\tau }^3+\frac{d}{6}{\tau }^3-\left(\frac{a^{2}b}{24}+\frac{\mathit{ad}}{8}\right){\tau }^4+\left(\frac{a^{3}b+7a^{2}d}{120}\right){\tau }^5-\left(\frac{a^{4}b}{720}+\frac{a^{3}d}{48}\right){\tau }^6$$ and

(4.12)

$$B(\tau )=-\tau +\frac{a}{2}{\tau }^2-\frac{a^2}{6}{\tau }^3+\frac{a^3}{24}{\tau }^4-\frac{a^4}{120}{\tau }^5+\frac{a^5}{720}{\tau }^6\text{.}$$

The polynomials $A(\tau )$ and $B(\tau )$ are the same as the Taylor expansion obtained above for $\alpha (\tau )$ and $\beta (\tau )$ , respectively. This means the limit of infinitely many terms (4.11) and (4.12) yields the exact solution (4.10). The accuracy of the scheme is measured using the following relative error

(4.13)

$$E=\frac{|\alpha (\tau )-A(\tau )|}{|\alpha (\tau )|}$$ where $\alpha (\tau )$ and $A(\tau )$ represent the exact and approximate solutions, respectively.

Table 1 illustrates the convergence of the HPTM. At $\tau =1$ and as j increases from 4 to 10 the error in $\alpha$ decreases from order ${10}^{-4}$ to ${10}^{-8}$ and from ${10}^{-2}$ to ${10}^{-11}$ for $\beta$ . At $\tau =0.1$ and as j increases from 4 to 10 the error in $\alpha$ decreases from order ${10}^{-5}$ to ${10}^{-14}$ and from ${10}^{-8}$ to ${10}^{-16}$ for $\beta$ .

Table 1 Convergence of HPTM of the Vasicek model for $a=0.5\text{,}b=0.3$ and $d=0.1$ .

		$\tau =0.1$	$\tau =0.2$	$\tau =0.5$	$\tau =0.7$	$\tau =1$
$j=4$	$A(\tau )$	1.19416E−5	9.66072E−5	1.55847E−3	4.36287E-3	1.30838E-2
	$B(\tau )$	5.29545E−8	8.61305E−7	3.53106E−5	1.39979E−4	6.10400E−4
$j=6$	$A(\tau )$	2.72792E−9	8.86011E−8	9.02825E−6	4.98717E−5	3.08174E−4
	$B(\tau )$	3.1584E−12	2.05922E−10	5.30832E−8	4.14078E−7	3.70637E−6
$j=8$	$A(\tau )$	3.84661E−13	4.88715E−11	3.13493E−8	3.40916E−7	4.32664E−6
	$B(\tau )$	8.13539E−17	2.92069E−14	4.63624E−11	7.1053E−10	1.30249E−8
$j=10$	$A(\tau )$	1.11102E−14	1.99155E−15	7.16094E−11	1.53097E−9	3.98281E−8
	$B(\tau )$	1.91192E−16	5.38435E−16	2.64043E−14	7.95826E−13	2.98415E−11

Fig. 1 (a) shows that the exact and numerical solutions of Eq. (4.1) are in good agreement. The bond price behaviour as a function of $\tau$ and $r(\tau )$ is captured in Fig. 1(b).

Example 4.2

In this experiment, we consider the Cox-Ingerson-Ross (CIR) model (see Björk, 2004). This model is obtained from Eq. (4.1) when the parameters are $a(\tau )=-a\text{,}b(\tau )=\mathit{ab}\text{,}c(\tau )={\sigma }^2$ and $d(\tau )=0$ . The exact solution of this model is

(4.14)

$$\alpha (\tau )=\frac{2\mathit{ab}}{{\sigma }^2}\ln \left[\frac{2\gamma e^{(\gamma +a)\frac{\tau }{2}}}{2\gamma +\left(1-e^{\gamma \tau }\right)\left(\gamma +a\right)}\right]\text{,}$$

(4.15)

$$\beta (\tau )=-\frac{2\left(e^{\gamma \tau }-1\right)}{2\gamma +\left(1-e^{\gamma \tau }\right)\left(\gamma +a\right)}\text{,}$$ where $\gamma =\sqrt{a{\sigma }^2+b^2}$ . The corresponding Taylor expansion is $$\begin{array}{c}\alpha (\tau )=-\frac{\mathit{ab}}{2}{\tau }^2+\frac{a^{2}b}{6}{\tau }^3+\left(\frac{\mathit{ab}{\sigma }^2}{24}-\frac{a^{3}b}{24}\right){\tau }^4+\left(\frac{a^{4}b}{120}-\frac{4a^{2}b{\sigma }^2}{120}\right){\tau }^5+O\left[{\tau }^6\right]\text{,}\hfill \\ \beta (\tau )=-\frac{a}{2}{\tau }^2+\left(\frac{a^2}{6}-\frac{{\sigma }^2}{6}\right){\tau }^3+\left(\frac{4a{\sigma }^2}{24}-\frac{a^3}{24}\right){\tau }^4+\left(\frac{a^4}{120}-\frac{11a^2{\sigma }^2}{120}+\frac{{\sigma }^4}{30}\right){\tau }^5+O\left[{\tau }^6\right]\text{.}\hfill \end{array}$$

Close

Fig. 1

(a) Exact and numerical solutions of 4.1 and (b) Bond price process for the Vasicek model when $a=0.5\text{,}b=0.3$ and $d=0.1$ at $j=8$ .

Export to PPT

Using the initial conditions ${\alpha }_0(\tau )=A(0)=\alpha (0)=0$ and ${\beta }_0(\tau )=B(0)=\beta (0)=u$ , we solve Eq. (4.1) for $A_j(\tau )\text{,}{B}_j(\tau )\text{,}j=0\text{,}1\text{,}\dots$ . Considering the bond that pay-off 1 at maturity implies $u=0$ . The solution computed using the HPTM at order $j=6$ is given by $$\begin{array}{c}\alpha (\tau )=-\frac{\mathit{ab}}{2}{\tau }^2+\frac{a^{2}b}{6}{\tau }^3+\left(\frac{\mathit{ab}{\sigma }^2}{24}-\frac{a^{3}b}{24}\right){\tau }^4+\left(\frac{a^{4}b}{120}-\frac{4a^{2}b{\sigma }^2}{120}\right){\tau }^5\text{,}\hfill \\ \beta (\tau )=-\frac{a}{2}{\tau }^2+\left(\frac{a^2}{6}-\frac{{\sigma }^2}{6}\right){\tau }^3+\left(\frac{4a{\sigma }^2}{24}-\frac{a^3}{24}\right){\tau }^4+\left(\frac{a^4}{120}-\frac{11a^2{\sigma }^2}{120}+\frac{{\sigma }^4}{30}\right){\tau }^5\text{.}\hfill \end{array}$$

Again we observe that the Taylor expansion and the solution of the CIR model computed by the HPTM are in good agreement. Table 2 records the error for different values of j and different values of $\tau$ taken randomly. The same conclusion applies as to the Vasicek model, that the error decreases rapidly as $\tau$ gets closer to 0. The bond price in terms of time to maturity $\tau$ and the interest rate r is given by

(4.16)

$$p(\tau \text{,}r)=e^{\alpha (\tau )-\beta (\tau )r(\tau )}\text{.}$$

Table 2 Convergence of the CIR model for $a=0.5\text{,}b=0.3\text{,}\sigma =0.1$ .

		$\tau =0.2$	$\tau =0.4$	$\tau =0.6$	$\tau =0.8$	$\tau =1$
$n=4$	$A(\tau )$	1.43085E−5	1.16921E−4	4.02833E−4	9.74221E−4	1.9403E−3
	$B(\tau )$	4.95614E−7	8.31425E−6	4.40686E−5	1.45623E−4	3.71224E−4
$n=6$	$A(\tau )$	1.07357E−10	5.13301E−9	5.16706E−8	2.71528E−7	9.93647E−7
	$B(\tau )$	2.00642E−10	1.30547E−8	1.51097E−7	8.62187E−7	3.33847E−6
$n=8$	$A(\tau )$	2.68126E−12	1.90102E−10	3.28522E−9	2.49695E−8	1.20746E−7
	$B(\tau )$	1.28341E−13	3.43621E−11	9.14227E−10	9.47091E−9	5.84918E−8
$n=10$	$A(\tau )$	1.22378E−12	7.96574E−13	5.11847E−12	7.03396E−11	5.42104E−10
	$B(\tau )$	7.29208 E−16	1.16398E−14	6.42225E−13	1.09941E−11	9.83115E−11

Fig. 2 shows the behaviour of the parameters $\alpha$ and $\beta$ as well as the corresponding bond price (4.16) for $a=0.5\text{,}b=0.3\text{,}\sigma =0.1$ . Clearly, exact and numerical solutions are in good agreement.

Close

Fig. 2

(a) Parameters $\alpha (\tau )$ and $\beta (\tau )$ exact and approximate and (b) bond price behaviour for $0⩽r⩽0.5\text{,}a=0.5\text{,}b=0.3\text{,}\sigma =0.1$ .

Export to PPT