Finite Difference Method for the Hull–White Partial Differential Equations

Abstract

This paper reviews the finite difference method (FDM) for pricing interest rate derivatives (IRDs) under the Hull–White Extended Vasicek model (HW model) and provides the MATLAB codes for it. Among the financial derivatives on various underlying assets, IRDs have the largest trading volume and the HW model is widely used for pricing them. We introduce general backgrounds of the HW model, its associated partial differential equations (PDEs), and FDM formulation for one- and two-asset problems. The two-asset problem is solved by the basic operator splitting method. For numerical tests, one- and two-asset bond options are considered. The computational results show close values to analytic solutions. We conclude with a brief comment on the research topics for the PDE approach to IRD pricing.

1. Introduction

The finite difference method (FDM) has been widely used to compute the prices of financial derivatives numerically (Duffy [1]). Many studies on option pricing employ FDM, especially on models such as Black–Scholes partial differential equation (PDE) for equity or exchange rate derivatives (Kim et al. [2]). Among financial derivatives on various underlying assets, interest rate derivatives (IRDs), whose payoffs depend on interest rates or bond prices, have the largest trading volume in the global over-the-counter market (BIS [3]). However, despite their popularity, FDM for IRD pricing has received relatively less attention.

The PDE approach to IRD pricing is naturally associated with short rate as interest rates are defined under a single probability measure and Markov in the models. The Hull–White extended Vasicek model (Hull and White [4]), also called Hull–White (HW) model, is the most widely used one-factor short rate model for pricing IRDs. Few studies have focused on finding numerical solutions under the HW model using FDM. For single-asset problems, Hull and White [5] suggested a modified explicit FDM under the HW model for better convergence. Vetzal [6] developed an implicit FDM to improve upon Hull and White [5]. In [7], Falco et al. proposed an algorithm for the complementarity problem arising from pricing American put options on a zero-coupon bond under the HW model by using $\theta$-method for discretization in time and LU decomposition. For multi-asset problems, Dang et al. [8] proposed a general PDE pricing framework for cross-currency IRDs under the two-asset HW model allowing a stochastic currency rate.

This paper reviews the HW model and the associated FDM formulation under a simplified setting for pricing IRDs. We also provide the MATLAB codes for one- and two-asset problems as Supplementary Materials. By doing so, we contribute to drawing the interest of researchers, especially beginners, to numerical PDE approaches to IRD pricing and saving their time to understand the topic.

The remainder of this paper is organized as follows. Section 2 reviews the single- and the two-asset HW models and PDEs. Section 3 formulates the procedure for solving the two-asset HW PDE using FDM. Section 4 presents numerical examples for one- and two-asset zero-coupon bond options. Section 5 concludes with comments.

2. Models and Partial Differential Equations

One advantage of the HW model that makes practitioners prefer it is analytical tractability. It provides the closed-form pricing formulas for bonds and standard IRDs such as caps/floors and swaptions. More generally, Turfus [9] derived a systematic method for exact analytic pricing formulas for European-style interest rate options under the HW model. Nevertheless, we cannot derive analytic pricing formulas under the HW model for many non-standard IRDs, especially with Bermudan option property, and thus, one should take a numerical approach. In this case, practitioners generally prefer the PDE approach to Monte Carlo simulation (MCS) for pricing derivatives since PDE returns the set of prices in the domain whereas MCS produces only one price point. This motivates the use of the HW PDEs and FDM to solve them.

This section briefly reviews interest rate models and introduces the HW PDEs for the single- and the two-asset IRDs. (Brigo and Mercurio [10] provide excellent reviews and in-depth analyses of interest rate models for pricing IRDs from both theoretical and practical perspectives).

2.1. Review of Interest Rate Models

Continuous-time interest rate models for pricing general (both plain vanilla and exotic) IRDs are of two major types. The first takes (hypothetical) instantaneous spot or forward rates under the risk-neutral measure as the basis for modeling the entire yield curve. Let $P(t,T)$ be the price at time t of the zero-coupon bond maturing at time T with the face value 1. The instantaneous spot rate, called “short rate”, is the interest rate $r\left(t\right)$ at time t over the infinitesimal time interval $[t,t+dt)$ used to calculate $P(t,T)$ by

$$P(t,T)=E_t\left[e^{-{\int }_t^{T}r\left(u\right)du}\right],$$

(1)

where $E_t[·]$ is the expectation calculated at time t under the risk-neutral measure. Short rate models define the stochastic dynamics of $r\left(t\right)$ (Black and Karasinski [11], Cox et al. [12], Hull and White [4], and Vasicek [13]). Once we have a short rate model, all interest rates on the yield curve can be derived as $r\left(t\right)$ dominates the entire yield curve. One of the advantages of the short rate models is that they often provide analytical formulas for the prices of bonds and liquidly traded IRDs (Wong and Zhao [14]). In addition, both the PDE approach and MCS can be directly applicable as $r\left(t\right)$ in short rate models is defined under a single probability measure and satisfies Markov property. However, they cannot produce a realistic volatility structure without introducing some complex assumptions. (One way for overcoming this drawback is two-factor modeling, for example, Hull and White [15]).

The instantaneous forward rate $f(t,T)$ at time t over the infinitesimal time interval $[T,T+dT)$ is defined to satisfy the following:

$$P(t,T)=E_t\left[e^{-{\int }_t^{T}f(u,T)du}\right],$$

(2)

where T is the maturity of $f(t,T)$. Instantaneous forward rate models, also known as Heath–Jarrow–Merton framework (Heath et al. [16]), define the stochastic dynamics of $f(t,T)$ for all T (i.e., the forward rate curve) under the single risk-neutral measure simultaneously (Heath et al. [16], Mercurio and Moraleda [17], and Ritchken and Sankarasubramanian [18]). They are theoretically more general than short rate models as they allow more flexible volatility specifications and take short rate models as a special case. However, they are usually non-Markov, and thus, we cannot derive the PDE for pricing IRDs directly. (The short rate from the Heath–Jarrow–Merton framework is Markov if the volatility function is deterministic and separable into a product of a function t and a function of T, as illustrated by Carverhill [19]).

Both the instantaneous interest rate models have two main common limitations: (1) they are expressed in terms of unobserved variables, (2) therefore, it is complicated to calibrate the model to the prices of caps/floors and swaptions. These shortcomings have led the second approach, namely the market model, which takes observable market interest rates, such as LIBOR or swap rates, as the basis for modeling and describes the evolution of their forward rates. If we model forward LIBORs (Brace et al. [20], Jamshidian [21], and Miltersen et al. [22]), we call the model “LIBOR market model” or “lognormal forward LIBOR model” or “Brace–Gatarek–Musiela model”. If we model forward swap rates (Jamshidian [21]), we call the model “swap market model” or “lognormal forward swap model”. This choice of variables for stochastic modeling in the market models makes the parameter calibration much more direct and simpler.

Market models assume that forward rates in the market follow lognormal distributions under their respective equivalent martingale measures. That is, the behaviors of all market rates are not modeled simultaneously under a single probability measure and different probability measures are used for each market rate. Therefore, MCS is usually preferred when pricing IRDs using the market models.

In conclusion, short rate models are the most suitable for the numerical PDE approach to general IRD pricing, which has the highest priority in the case that closed-form pricing formulas are not available, and thus, remains popular even after more advanced models are developed. The HW model is the most widely used one-factor short rate model because of its no-arbitrage property and analytical tractability. In the next subsections, we introduce one- and two-asset HW models and their associated PDEs.

2.2. HW PDE for a Single Interest Rate

Short rate models are classified into two categories. One is “equilibrium model” (Cox et al. [12] and Vasicek [13]) which starts with the assumptions on economic variables and derives the behavior of interest rates. In the equilibrium model, the yield curve observed in the market does not coincide with the yield curve derived from the model since it is not an input but an output of the model. The other is “no-arbitrage model”. “No-arbitrage” means that the model parameters are consistent with the bond prices or interest rates observed in the market. The HW model is the most widely used one-factor short rate model because of its no-arbitrage property and analytical tractability.

Hull and White [4] extended the Vasicek’s [13] model so that it has time-varying parameters fitting the yield curve observed in the market to satisfy the no-arbitrage condition. In the HW model, the stochastic dynamics of a short rate $r\left(t\right)$ is assumed to be the following Ornstein–Uhlenbeck process:

$$dr\left(t\right)=\left[\theta \left(t\right)-ar\left(t\right)\right]dt+\sigma dW\left(t\right)$$

(3)

with $r\left(0\right)=r_0$, where $W\left(t\right)$ is the standard Brownian Motion under the risk-neutral measure. The function $\theta \left(t\right)=\frac{\partial f^M(0,t)}{\partial t}+a{f}^M(0,t)+\frac{{\sigma }^2}{2a}(1-e^{-2at})$ is chosen so that the model exactly replicates the current yield curve, where $f^M(0,T)=-\frac{\partial P^M(0,T)}{\partial T}$ and $P^M(0,T)$ is the current market discount factor for the maturity T. The parameters $\sigma$ and a define the volatility of the short rate in the future and the relative volatilities of long-term and short-term interest rates, respectively. In fact, Hull and White [4] also included time-dependent volatility parameters $a\left(t\right)$ and $\sigma \left(t\right)$ in the model to fit the option prices observed in the market. However, Hull and White [23] discussed that the time-dependent volatility parameters could exhibit non-stationarity and lead to unacceptable results when used to price exotic instruments. Consequently, Hull and White [24] recommend keeping the volatility parameters constant, which is Equation (3).

The extant literature uses two types of HW PDEs. One is the original PDE with respect to $\left(r\right(t),t)$, generally presented in the literature (Dang et al. [8], Hull and White [4,5], and Vetzal [6]). The other takes $x\left(t\right)$, the random component of $r\left(t\right)$, as the state variable from the following decomposition:

$$r\left(t\right)=\alpha \left(t\right)+x\left(t\right)$$

(4)

with $x\left(0\right)=0$, where $d\alpha \left(t\right)=\left[\theta \left(t\right)-a\alpha \left(t\right)\right]dt$ and $dx\left(t\right)=-ax\left(t\right)dt+\sigma dW\left(t\right)$ (Falco et al. [7] and Wong and Zhao [14]). Notice that $x\left(t\right)$ governs the randomness of $r\left(t\right)$ and $\alpha \left(t\right)$ is the deterministic component to make the HW model fit the discount factors (or zero-coupon rates) implied in the market yield curve. This approach is motivated by Hull and White’s [25] two-step trinomial tree model, which is an improvement upon Hull and White [26]. We take the latter as Falco et al. [7] demonstrated that it improves computation times and accuracy across various discretization methods. (Our formulation is not exactly the same as Falco et al. [7]).

Let $u(x,t)$ be the value at time t of a contingent claim on a single interest rate with $u(x,T)=\mathsf{\Phi }\left(x\right)$. We consider the following HW PDE with respect to $x\left(t\right)$ in Equation (4) and t:

$$\frac{\partial u(x,t)}{\partial t}=ax\frac{\partial u(x,t)}{\partial x}-\frac{1}{2}{\sigma }^2\frac{{\partial }^{2}u(x,t)}{\partial x^2}+[\alpha \left(t\right)+x]u(x,t)$$

(5)

for $(x,t)\in \mathbb{R}\times [0,T)$, where $\alpha \left(t\right)=f^M(0,t)+\frac{{\sigma }^2}{2a^2}{(1-e^{-at})}^2$. Equation (5) can be derived by the usual no-arbitrage argument applying Ito’s lemma.

2.3. HW Model and PDE for Two Interest Rates

Multi-asset IRDs are actively traded in the market as a form of non-standard swaps. The typical case is the quanto swap, whose underlying assets are in different currencies with payments made in the same currency. The nominal amount is usually denominated in one party’s domestic currency. For example, “differential swaps” (Chang et al. [27]) exchange two cash flows, one of which is determined by the domestic interest rate and the other is determined by the foreign interest rate. A “dual range accrual swap” is a quanto swap in which one of the parties pays the daily accrued interest determined by both the domestic and the foreign interest rates. If both underlying interest rates stay in their respective ranges on a date, the daily interest is accrued.

Table 1. An example of KRW dual range accrual swap.

Condition	Detail
Nominal amount	KRW 100,000,000
Maturity	15 years
Payment currency	KRW
Interest 1	$5\%\times \frac{n}{365}$, Annually
	n = number of calendar days satisfying the following conditions:
	0 < KRW 91-day CD rate < 5% and 0 < USD 3-month LIBOR < 5%
Interest 2	KRW 91-day CD rate + 0.05%, Quarterly
Cancellable option	The interest 1 payer has the right to cancel the contract every year from trade.

exhibits an example of KRW dual range accrual swaps whose nominal amount is denominated in KRW and all payments are made in KRW, but the underlying interest rates include USD 3-month LIBOR.

Closed-form pricing formulas for some quanto swaps could be available under the HW model even if they are non-standard swaps. For example, Brigo and Mercurio [10] found the closed-form formulas for differential swaps under the G2++ model, which takes the HW model as a special case. However, a large amount of non-standard IRDs traded in the market involve Bermudan optionality like the cancellable option in Table 1, and hence, we should take a numerical approach for pricing. This motivates our use of the two-asset HW PDE and FDM to solve it.

The HW model and PDE for contingent claims on a single interest rate in Section 2.2 can be extended to multi-asset problems. We consider the IRDs on two-currency assets, which consist of domestic and foreign interest rates. The payment is assumed to be made in the domestic currency, similar to the example in Table 1. We suppose that each currency has only one yield curve following the HW model. This means that the yield curves for calculating future cash flows and discounting are the same (The framework using separate yield curves for future cash flows and discounting for pricing IRDs with a single currency is called “multi-curve framework”. See Henrard [28] for more details). For the exchange rate between the two currencies, We consider the Geometric Brownian Motion satisfying the no-arbitrage condition between the two economies (Dang et al. [8] and Piterbarg [29]).

Let $r_1\left(t\right)$ and $r_2\left(t\right)$ be the short rates for the domestic and the foreign yield curves, respectively. Then the evolutions of $r_1\left(t\right)$ and $r_2\left(t\right)$ under the domestic risk-neutral measure are

$$\begin{array}{ccc}{r}_1\left(t\right)& =& {\alpha }_1\left(t\right)+x\left(t\right)\hfill \\ r_2\left(t\right)& =& {\alpha }_2\left(t\right)+y\left(t\right)\hfill \\ dx\left(t\right)& =& -a_{1}x\left(t\right)dt+{\sigma }_{1}d{W}_1\left(t\right)\hfill \\ dy\left(t\right)& =& -[{\rho }_{23}{\sigma }_2{\sigma }_3+a_{2}y\left(t\right)]dt+{\sigma }_{2}d{W}_2\left(t\right)\hfill \\ {\displaystyle \frac{dS\left(t\right)}{S\left(t\right)}}& =& [r_1\left(t\right)-r_2\left(t\right)]dt+{\sigma }_{3}d{W}_3\left(t\right)\hfill \end{array}$$

(6)

with ${\alpha }_i\left(t\right)=f_i^M(0,t)+\frac{{\sigma }_i^2}{2a_i^2}{(1-e^{-a_{i}t})}^2$ and $f_i^M(0,T)=-\frac{\partial P_i^M(0,T)}{\partial T}$ for $i=1,2$, where $P_i^M(0,T)$ is the price of the zero-coupon bond maturing at time T implied in the corresponding yield curve. The standard Brownian Motions $W_i\left(t\right)$, $W_2\left(t\right)$, and $W_3\left(t\right)$ are defined under the domestic risk-neutral measure with $d{W}_i\left(t\right)d{W}_j\left(t\right)={\rho }_{ij}dt$ for $i=1,2,3$. (The constant term $-{\rho }_{23}{\sigma }_2{\sigma }_3$ in the drift of $dy\left(t\right)$ is the adjustment to describe the foreign interest rate $y\left(t\right)$ under the domestic risk-neural measure. This is called “quanto adjustment”. If we set ${\rho }_{23}=0$ or ${\sigma }_3=0$, we can ignore this adjustment). Both $r_1\left(t\right)$ and $r_2\left(t\right)$ are used for calculating the future prices of underlying assets. However, as the payoff is settled in the domestic currency, $r_1\left(t\right)$ is applied to discounting.

Let $u(x,y,t)$ be the value of a contingent claim on the two interest rates at time t satisfying $u(x,y,T)=\mathsf{\Phi }(x,y)$. We can then derive the following two-asset HW PDE with respect to $x\left(t\right)$, $y\left(t\right)$ in Equation (6) and t:

$$\begin{array}{ccc}& {\displaystyle \frac{\partial u(x,y,t)}{\partial t}}=& a_{1}x{\displaystyle \frac{\partial u(x,y,t)}{\partial x}}-{\displaystyle \frac{{\sigma }_1^2}{2}}{\displaystyle \frac{{\partial }^{2}u(x,y,t)}{\partial x^2}}\hfill \\ & & +({\rho }_{23}{\sigma }_2{\sigma }_3+a_{2}y){\displaystyle \frac{\partial u(x,y,t)}{\partial y}}-{\displaystyle \frac{{\sigma }_2^2}{2}}{\displaystyle \frac{{\partial }^{2}u(x,y,t)}{\partial y^2}}\hfill \\ & & -{\rho }_{12}{\sigma }_1{\sigma }_2{\displaystyle \frac{{\partial }^{2}u(x,y,t)}{\partial x\partial y}}+[{\alpha }_1\left(t\right)+x]u(x,y,t)\hfill \end{array}$$

(7)

for $(x,y,t)\in {\mathbb{R}}^2\times [0,T)$. Please note that the coefficient function $[{\alpha }_1\left(t\right)+x]$ of the last term in the right-hand side is equal to $r_1$ at time t in the discounting curve.

3. Finite Difference Method

In this section, we formulate the FDM to solve the two-asset HW PDE in Equation (7) applying operator splitting method (OSM). The single-asset HW PDE in Equation (5) can be solved similarly.

The choice of a finite difference scheme is important when solving a PDE numerically. Dang et al. [8] numerically found that the Hundsdorfer and Verwer Alternating Direction Implicit method is more efficient than the Crank-Nicolson method armed with the Generalized Minimal Residual method, whose preconditioner is solved by the Fast Fourier Transform for the two-asset HW PDE with a stochastic exchange rate. (To the best of our knowledge, Dang et al. [8] is the only published paper on pricing multi-asset IRDs by the HW PDE). Jeong and Kim [30] provided computational results illustrating that the basic OSM provides better accuracy and robustness than the basic Alternating Direction Implicit method, although this was for two-asset option pricing problems in the Black-Scholes model and not for the two-asset HW PDE case. Based on these findings, we apply the basic OSM similar to Jeong and Kim [30]. Please see Hout and Toivanen [31] for more advanced OSMs.

We define a spatial differential operator $\mathcal{L}$ as

$$\mathcal{L}u=a_{1}x{\displaystyle \frac{\partial u}{\partial x}}+[{\rho }_{23}{\sigma }_2{\sigma }_3+a_{2}y]{\displaystyle \frac{\partial u}{\partial y}}-{\displaystyle \frac{1}{2}}{\sigma }_1^2{\displaystyle \frac{{\partial }^{2}u}{\partial x^2}}-{\displaystyle \frac{1}{2}}{\sigma }_2^2{\displaystyle \frac{{\partial }^{2}u}{\partial y^2}}-{\rho }_{12}{\sigma }_1{\sigma }_2{\displaystyle \frac{{\partial }^{2}u}{\partial x\partial y}}+[{\alpha }_1\left(t\right)+x]u.$$

(8)

Then the HW PDE in Equation (7) can be written as

$$\frac{\partial u}{\partial t}=\mathcal{L}u$$

(9)

for $(x,y,t)\in {\mathbb{R}}^2\times [0,T)$ with the terminal condition $u(x,y,T)=\mathsf{\Phi }(x,y)$. To apply the OSM, we define the following fractional differential operators:

$$\begin{array}{ccc}\hfill {\mathcal{L}}_{x}u& =& a_{1}x{\displaystyle \frac{\partial u}{\partial x}}-{\displaystyle \frac{1}{2}}{\sigma }_1^2{\displaystyle \frac{{\partial }^{2}u}{\partial x^2}}-{\displaystyle \frac{1}{2}}{\rho }_{12}{\sigma }_1{\sigma }_2{\displaystyle \frac{{\partial }^{2}u}{\partial x\partial y}}+{\displaystyle \frac{1}{2}}[{\alpha }_1\left(t\right)+x]u\hfill \\ \hfill {\mathcal{L}}_{y}u& =& [{\rho }_{23}{\sigma }_2{\sigma }_3+a_{2}y]{\displaystyle \frac{\partial u}{\partial y}}-{\displaystyle \frac{1}{2}}{\sigma }_2^2{\displaystyle \frac{{\partial }^{2}u}{\partial y^2}}-{\displaystyle \frac{1}{2}}{\rho }_{12}{\sigma }_1{\sigma }_2{\displaystyle \frac{{\partial }^{2}u}{\partial x\partial y}}+{\displaystyle \frac{1}{2}}[{\alpha }_1\left(t\right)+x]u.\hfill \end{array}$$

(10)

Then $\mathcal{L}$ is split into $\mathcal{L}={\mathcal{L}}_x+{\mathcal{L}}_y$ and we can obtain a numerical solution for Equation (9) by using this splitting scheme.

To solve Equation (9) numerically, we truncate its unbounded domain into bounded one as $(x,y,t)\in (x_{min},x_{max})\times (y_{min},y_{max})\times [0,T)$ for sufficiently small $x_{min}$, $y_{min}$ and large $x_{max}$, $y_{max}$ so that the region outside the truncated domain hardly affects our numerical results. We then discretize the finite domain uniformly as $\left(x_i,y_j,t_k\right)=\left(x_{min}+(i-1)h_x,y_{min}+(j-1)h_y,(k-1)h_t\right)$ for $i=1,\cdots ,N_x+1$, $j=1,\cdots ,N_y+1$, and $k=1,\cdots ,N_t+1$, where $h_x=\frac{x_{max}-x_{min}}{N_x}$, $h_y=\frac{y_{max}-y_{min}}{N_y}$ and $h_t=\frac{T}{N_t}$. Please note that we start our indices not from 0 but 1 to match them with MATLAB’s index rule.

Let $u_{i,j}^k$ be the approximated value of $u(x_i,y_j,t_k)$ from FDM. By applying central difference in space and implicit difference in time, ignoring error, we have the following splitting scheme:

$$\begin{array}{cc}\hfill \frac{u_{i,j}^{k+1}-u_{i,j}^{k+\frac{1}{2}}}{h_t}& ={\left({\mathcal{L}}_{x}u\right)}_{i,j}^{k+\frac{1}{2}}\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill \frac{u_{i,j}^{k+\frac{1}{2}}-u_{i,j}^k}{h_t}& ={\left({\mathcal{L}}_{y}u\right)}_{i,j}^k\hfill \end{array}$$

(12)

for $i=1,\cdots ,N_x-1$, $j=1,\cdots ,N_y-1$, and $k=0,\cdots ,N_t-1$. Here, the discretized fractional operators are defined as

$$\begin{array}{cc}\hfill {\left({\mathcal{L}}_{x}u\right)}_{i,j}^{k+\frac{1}{2}}& =a_1{x}_i{D}_x{u}_{i,j}^{k+\frac{1}{2}}-\frac{{\sigma }_1^2}{2}{D}_{xx}{u}_{i,j}^{k+\frac{1}{2}}-\frac{{\rho }_{12}{\sigma }_1{\sigma }_2}{2}{D}_{xy}{u}_{i,j}^{k+1}+\frac{1}{2}[{\alpha }_1\left(t_{k+1}\right)+x_i]u_{i,j}^{k+1}\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill {\left({\mathcal{L}}_{y}u\right)}_{i,j}^k& =[{\rho }_{23}{\sigma }_2{\sigma }_3+a_2{y}_j]D_y{u}_{i,j}^k-\frac{{\sigma }_2^2}{2}{D}_{yy}{u}_{i,j}^k-\frac{{\rho }_{12}{\sigma }_1{\sigma }_2}{2}{D}_{xy}{u}_{i,j}^{k+\frac{1}{2}}+\frac{1}{2}[{\alpha }_1(t_{k+\frac{1}{2}})+x_i]u_{i,j}^{k+\frac{1}{2}}\hfill \end{array}$$

(14)

for a fixed time index k, where the first and the second order differences are defined as

$$\begin{array}{cc}\hfill D_x{u}_{i,j}& =\frac{u_{i+1,j}-u_{i-1,j}}{2h_x},D_{xx}{u}_{i,j}=\frac{u_{i+1,j}+2u_{i,j}-u_{i-1,j}}{h_x^2}\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill D_y{u}_{i,j}& =\frac{u_{i,j+1}-u_{i,j-1}}{2h_y},D_{yy}{u}_{i,j}=\frac{u_{i,j+1}+2u_{i,j}-u_{i,j-1}}{h_y^2}\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill D_{xy}{u}_{i,j}& =\frac{u_{i+1,j+1}-u_{i+1,j-1}-u_{i-1,j+1}+u_{i-1,j-1}}{4h_x{h}_y}\hfill \end{array}$$

(17)

for $i=1,\cdots ,N_x-1$ and $j=1,\cdots ,N_y-1$. We use the linear boundary conditions $\frac{{\partial }^{2}u(x_{min},y,t)}{\partial x^2}=\frac{{\partial }^{2}u(x_{max},y,t)}{\partial x^2}=\frac{{\partial }^{2}u(x,y_{min},t)}{\partial y^2}=\frac{{\partial }^{2}u(x,y_{max},t)}{\partial y^2}=0$ to be discretized into $u_{1,j}=2u_{2,j}-u_{3,j}$, $u_{N_x+1j}=-u_{N_x-1,j}+2u_{N_x,j}$, $u_{i,1}=2u_{i,2}-u_{i,3}$, and $u_{i,N_y+1}=-u_{i,N_y-1}+2u_{i,N_y}$. Please note that Equations (11) and (13) go from the known values at time $t_{k+1}$ to unknown values at time $t_{k+\frac{1}{2}}$ and Equations (12) and (14) go from the known values at time $t_{k+\frac{1}{2}}$ to unknown values at time $t_k$.

We can derive a system of linear equations to solve Equation (11) from Equations (13), (15), and (17) along with the linear boundary conditions as

$$\begin{array}{ccccc}& \left[\begin{array}{cccccc}e+2d_2& f_2-d_2& 0& \cdots & 0& 0\\ d_3& e& f_3& \cdots & 0& 0\\ 0& d_4& e& \cdots & 0& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& 0& d_{N_x-1}& e& f_{N_x-1}\\ 0& 0& 0& \cdots & d_{N_x}-f_{N_x}& e+2f_{N_x}\end{array}\right]& \left[\begin{array}{c}{u}_{2,j}^{k+\frac{1}{2}}\\ u_{3,j}^{k+\frac{1}{2}}\\ u_{4,j}^{k+\frac{1}{2}}\\ ⋮\\ u_{N_x-1,j}^{k+\frac{1}{2}}\\ u_{N_x,j}^{k+\frac{1}{2}}\end{array}\right]& =& \left[\begin{array}{c}{g}_{2,j}^{k+1}\\ g_{3,j}^{k+1}\\ g_{4,j}^{k+1}\\ ⋮\\ g_{N_x-1,j}^{k+1}\\ g_{N_x,j}^{k+1}\end{array}\right]\\ \\ (& \mathbf{A}& {\mathbf{u}}_{2:N_x,j}^{k+\frac{1}{2}}& =& {\mathbf{g}}_{2:N_x,j}^{k+1}& )\end{array}$$

(18)

for fixed indices j and k. Here, the elements of the coefficient matrix $\mathbf{A}$ in the linear system (18) are defined as

$$\begin{array}{c}{d}_i=-p_i-q,e=1+2q,f_i=p_i-q,\\ g_{i,j}^{k+1}=s\left(u_{i+1,j+1}^{k+1}-u_{i+1,j-1}^{k+1}-u_{i-1,j+1}^{k+1}+u_{i-1,j-1}^{k+1}\right)+\left(1-w^{k+1}\right)u_{i,j}^{k+1},\end{array}$$

(19)

where $p_i=a_1{\displaystyle \frac{h_t}{2h_x}}{x}_i$, $q={\displaystyle \frac{1}{2}}{\sigma }_1^2{\displaystyle \frac{h_t}{h_x^2}}$, $s={\displaystyle \frac{1}{2}}{\rho }_{12}{\sigma }_1{\sigma }_2{\displaystyle \frac{h_t}{4h_x{h}_y}}$, and $w_i^{k+1}={\displaystyle \frac{1}{2}}\left[x_i+\alpha \left(t_{k+1}\right)\right]h_t$.

Similarly, Equation (12) to find $u_{i,j}^k$s from $u_{i,j}^{k+\frac{1}{2}}$s also can be expressed in a system of linear equations as

$$\begin{array}{ccccc}& \left[\begin{array}{cccccc}e+2d_2& f_2-d_2& 0& \cdots & 0& 0\\ d_3& e& f_3& \cdots & 0& 0\\ 0& d_4& e& \cdots & 0& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& 0& d_{N_y-1}& e& f_{N_y-1}\\ 0& 0& 0& \cdots & d_{N_y}-f_{N_y}& e+2f_{N_y}\end{array}\right]& \left[\begin{array}{c}{u}_{i,2}^k\\ u_{i,3}^k\\ u_{i,4}^k\\ ⋮\\ u_{i,N_y-1}^k\\ u_{i,N_y}^k\end{array}\right]& =& \left[\begin{array}{c}{g}_{i,2}^{k+\frac{1}{2}}\\ g_{i,3}^{k+\frac{1}{2}}\\ g_{i,4}^{k+\frac{1}{2}}\\ ⋮\\ g_{i,N_y-1}^{k+\frac{1}{2}}\\ g_{i,N_y}^{k+\frac{1}{2}}\end{array}\right]\\ \\ (& \mathbf{B}& {\mathbf{u}}_{i,2:N_y}^k& =& {\mathbf{g}}_{i,2:N_y}^{k+\frac{1}{2}}& )\end{array}$$

(20)

for fixed indices i and k. Here, the elements of the coefficient matrix $\mathbf{B}$ in the linear system (20) are defined as

$$\begin{array}{c}{d}_j=-p_j-q,e=1+2q,f_j=p_j-q,\\ g_{i,j}^{k+\frac{1}{2}}=s\left(u_{i+1,j+1}^{k+\frac{1}{2}}-u_{i+1,j-1}^{k+\frac{1}{2}}-u_{i-1,j+1}^{k+\frac{1}{2}}+u_{i-1,j-1}^{k+\frac{1}{2}}\right)+\left(1-w_i^{k+\frac{1}{2}}\right)u_{i,j}^{k+\frac{1}{2}},\end{array}$$

(21)

where $p_j=[{\rho }_{23}{\sigma }_2{\sigma }_3+a_2{y}_j]{\displaystyle \frac{h_t}{2h_y}}$, $q={\displaystyle \frac{1}{2}}{\sigma }_2^2{\displaystyle \frac{h_t}{h_y^2}}$, $s={\displaystyle \frac{1}{2}}{\rho }_{12}{\sigma }_1{\sigma }_2{\displaystyle \frac{h_t}{4h_x{h}_y}}$, and $w_i^{k+\frac{1}{2}}={\displaystyle \frac{1}{2}}\left[x_i+\alpha (t_{k+\frac{1}{2}})\right]h_t$. Please note that $w_i^{k+\frac{1}{2}}$ is not indexed by j.

4. Numerical Experiments

We consider numerical experiments for single- and two-asset at-the-money zero-coupon bond options. We assume that the domestic and the foreign zero-coupon rate curves are given as

Table 2. Zero rate curves for numerical tests. All rates are continuously compounded with actual/365 yield convention.

Time		Curve 1 (Domestic)		Curve 2 (Foreign)
Maturity	Days	Rate (%)	${\mathit{P}}_{\mathbf{1}}^{\mathit{M}}(\mathbf{0},\mathit{T})$	Rate (%)	${\mathit{P}}_{\mathbf{2}}^{\mathit{M}}(\mathbf{0},\mathit{T})$
1 day	1	3.03987	0.999917	0.12826	0.999996
3 months	93	3.44064	0.991272	0.25260	0.999357
6 months	184	3.50462	0.982488	0.40222	0.997974
9 months	275	3.54022	0.973680	0.32367	0.997564
1 year	369	3.58056	0.964449	0.36307	0.996336
2 years	733	3.66949	0.928958	0.60594	0.987905
3 years	1097	3.75130	0.893379	0.99992	0.970395
4 years	1462	3.80810	0.858531	1.43411	0.944176
5 years	1828	3.86056	0.824197	1.85815	0.911139
7 years	2560	3.95699	0.757652	2.56285	0.835479
10 years	3654	4.07967	0.664704	3.24935	0.722317

. Let $r_1$ and $r_2$ denote the short rates for the domestic and the foreign yield curves, respectively. For model parameters, we use $a_1=0.02$, $a_2=0.04$, ${\sigma }_1=0.008$, ${\sigma }_2=0.012$, ${\rho }_{12}=0.6$, and ${\rho }_{23}=0$. Regarding the parameters for FDM, we set $x_{min}=y_{min}=-0.20$, $x_{max}=y_{max}=0.20$, and $N_x=N_y=300$. We choose $N_t$ so that the grid size $h_t$ is equal to two calendar days since we will test various maturities.

The errors of our numerical solutions are measured by the root mean squared errors (RMSE) defined as

$$\begin{array}{cc}\hfill RMSE& =\sqrt{{\displaystyle \frac{1}{N}}\sum _i\sum _j{\left(u_{i,j}^0-u(x_i,y_j,0)\right)}^2}\hfill \end{array}$$

(22)

for $(x_i,y_j)\in \Omega$, where $\Omega =({\displaystyle \frac{x_{min}}{4}},{\displaystyle \frac{x_{max}}{4}})\times ({\displaystyle \frac{y_{min}}{4}},{\displaystyle \frac{y_{max}}{4}})$ and N is the total number of grid points in $\Omega$. We measure the errors in the sub-region $\Omega$ as $(x_i,y_j)$ in the area produces zero-coupon bond prices near the strike prices and large differences between the approximated and exact closed-form prices.

4.1. Single-Asset Options

Consider a domestic zero-coupon bond maturing at time S with a face value 1. The payoff function of the European call option with strike price K and maturity T ($T<S$) on this bond is

$$V\left(P(T,S)\right)={\left[P(T,S)-K\right]}^{+}.$$

(23)

If we calculate Equation (1) in the HW model, $P(t,S)$ is obtained as (Brigo and Mercurio [10])

$$P(t,S)=A(t,S)e^{-B(t,S)\left[\alpha \left(t\right)+x\left(t\right)\right]},$$

(24)

where $A(t,S)={\displaystyle \frac{P^M(0,S)}{P^M(0,t)}}{e}^{B(t,S)f^M(0,t)-\frac{{\sigma }^2}{4a}(1-e^{-2at})B{(t,S)}^2}$ and $B(t,S)={\displaystyle \frac{1}{a}}\left[1-e^{-a(S-t)}\right]$. Since $P(t,S)$ is a function of $x\left(t\right)$, Equation (23) can be rewritten as

$$\mathsf{\Phi }\left(x\right)={\left[A(T,S)e^{-B(T,S)\left[\alpha \left(T\right)+x\right]}-K\right]}^{+},$$

(25)

where $x=x\left(T\right)$. This expression is necessary to solve the HW PDE with a terminal condition with respect to x.

The closed-form formula for the European vanilla call option with the payoff in Equation (23) under the HW model is (Brigo and Mercurio [10])

$$\begin{array}{ccc}\hfill u(x,t)& =& ZBC(x,t,T,S,K)\hfill \\ \hfill & =& P(t,S)N\left(h\right)-KP(t,T)N(h-{\sigma }_P),\hfill \end{array}$$

(26)

where $x=x\left(t\right)$ with $x\left(0\right)=0$, ${\sigma }_P=\sigma \sqrt{{\displaystyle \frac{1-e^{-2a(T-t)}}{2a}}}B(T,S)$, $h={\displaystyle \frac{1}{{\sigma }_P}}\ln {\displaystyle \frac{P(t,S)}{P(t,T)K}}+{\displaystyle \frac{{\sigma }_P}{2}}$, and $N(·)$ is the cumulative distribution function of $N(0,1)$.

We choose T from 1–7 years in Table 2, $S=T+2$, and $K=P(0,S)$ which is the at-the-money strike.

Table 3. Computational results of the one-asset numerical experiments.

T (Year)	S (Year)	K	Numerical Solution	RMSE
1	3	0.893186925	0.031802834	0.0000088634
2	4	0.858436372	0.060984307	0.0000128773
3	5	0.824290518	0.087883163	0.0000145132
4	6	0.790788386	0.111866668	0.0000142805
5	7	0.757829933	0.133221129	0.0000121528
7	9	0.694856537	0.168386157	0.0000003208

illustrates the results of the numerical computations for the various T, S, and K.

4.2. Multi-Asset Options

For the multi-asset example, we consider a two-asset digital call option with strike prices $K_1$ and $K_2$ maturing at time T on two zero-coupon bonds with face values 1. Its payoff at maturity is expressed as

$$V\left(P_1(T,S_1),P_2(T,S_2)\right)=1_{\left\{P_1(T,S_1)\ge K_1\mathrm{and}{P}_2(T,S_2)\ge K_2\right\}}.$$

(27)

Since $P_1(T,S_1)$ and $P_2(T,S_1)$ are functions of $x\left(t\right)$ and $y\left(t\right)$, respectively, Equation (27) can be rewritten as

$$\mathsf{\Phi }(x,y)=1_{\left\{A_1(T,S_1)e^{-B_1(T,S_1)\left[{\alpha }_1\left(t\right)+x\right]}\ge K_1\mathrm{and}{A}_2(T,S_2)e^{-B_2(T,S_2)\left[{\alpha }_2\left(t\right)+y\right]}\ge K_2\right\}},$$

(28)

where $x=x\left(T\right)$ and $y=y\left(T\right)$. This expression is necessary to solve the two-asset HW PDE with a terminal condition with respect to x and y.

The closed-form formula for the two-asset digital call option in Equation (27) can be derived under the HW model as follows:

$$\begin{array}{ccc}\hfill u(x,y,t)& =& ZBDC(x,y,t,T,S_1,S_2,K_1,K_2)\hfill \\ & =& P_1(t,T)M(h_1-{\sigma }_{P_1},h_2-{\sigma }_{P_2};{\rho }_{12}),\hfill \end{array}$$

where $x=x\left(t\right)$ and $y=y\left(t\right)$ with $x\left(0\right)=y\left(0\right)=0$, ${\sigma }_{P_i}={\sigma }_i\sqrt{{\displaystyle \frac{1-e^{-2a_i(T-t)}}{2a_i}}}{B}_i(T,S_i)$, $h_i={\displaystyle \frac{1}{{\sigma }_{P_i}}}\ln {\displaystyle \frac{P_i(t,S_i)}{P_i(t,T)K_i}}+{\displaystyle \frac{{\sigma }_{P_i}}{2}}$ for $i=1,2$, and $M(·,·;\rho )$ is the cumulative distribution function of the bivariate standard normal distribution with correlation $\rho$.

Similar to the one-asset examples, we choose T from 1–7 years in Table 2, $S_1=S_2=T+2$, and $K_i=P_i(0,S_i)$ for $i=1,2$.

Table 4. Computational results of the two-asset numerical experiments.

T (Year)	${\mathit{S}}_1$ (Year)	${\mathit{S}}_2$ (Year)	${\mathit{K}}_1$	${\mathit{K}}_2$	Numerical Solution	RMSE
1	3	3	0.893186925	0.970272123	0.547677056	0.000388024
2	4	4	0.858436372	0.944094608	0.600037268	0.000459376
3	5	5	0.824290518	0.911237792	0.702678318	0.000500082
4	6	6	0.790788386	0.875784156	0.776749239	0.000517691
5	7	7	0.757829933	0.835708511	0.801148847	0.000474142
7	9	9	0.694856537	0.761627292	0.757202620	0.000321140

illustrates the results of the numerical computations for the various T, $S_i$s, and $K_i$s.

5. Conclusions

This paper reviews the HW model, its PDE, and the associated FDM for finding numerical prices of IRDs along with the related literature. In the FDM, the HW PDE is discretized uniformly by central differencing and implicitly in time. In the numerical examples, we compute prices for one- and two-asset zero-coupon bond options. The two-asset HW PDE is solved by OSM.

As the PDE approach to IRD pricing is relatively scant in the extant literature, many research questions on the HW PDE have hitherto not been investigated sufficiently. For example, the HW PDE for barrier options involves curved boundary problems as interest rates and bond prices are non-linear functions of short rates. In addition, the range accrual products shown in Table 1 have series of singular points for the same reason. Therefore, it is better to apply some remedy to spatial grids for obtaining stable numerical prices. In terms of a time-grid, IRDs usually have very long-term maturities. Hence, FDM for them needs an efficient time-grid structure to expedite the computation time. Lastly, in this paper, we used the basic OSM for review. A more efficient (higher-order) scheme can also be considered.