Unbiased estimators for the Heston model with stochastic interest rates

The classical Heston model (Heston [13]) is one of the fundamental models in finance, and it has been widely applied in various financial markets, such as the equity, fixed income and foreign exchange markets, due to its tractability in modelling the term structure of implied volatility. However, in this model, the interest rate is constant, and this assumption is often not appropriate for long-maturity options, because the long-term behaviour of the interest rate is typically far from constant. This phenomenon has been empirically investigated by Bakshi, Cao and Chen [2]. A natural extension of the Heston model is to incorporate a stochastic interest rate, which is usually referred to as the Heston model with stochastic interest rates. There are several contributions in this direction, such as Grzelak and Oosterlee [12], Van Haastrecht and Pelsser [21] and references therein.

Under the Heston model with stochastic interest rates, the price of an option can be written as

where $S$ is the solution to the Heston model with stochastic interest rates and $r_{t}$ is the underlying interest rate. Here, $P$ denotes the payoff functional of an option. We are interested in calculating this expectation, as for the majority of options, there are no closed-form formulas. A classical approach is to use a Monte Carlo method associated with a time-discrete scheme on $S$ for an approximate value. Recently, Rhee and Glynn [20] proposed several unbiased Monte Carlo estimators based on a randomization idea for a stochastic differential equation, which are unbiased versions of multilevel Monte Carlo estimators by Giles [9]. These unbiased estimators have a clear advantage over the standard Monte Carlo estimator, as the latter is typically biased when $S$ has to be approximated through a time-discrete scheme. To combine Rhee and Glynn’s estimators and the Heston model with stochastic interest rates, it is essential to develop a numerical scheme with a good convergence rate in the $L^{2}$ norm. However, standard theorems, such as those in Kloeden and Platen [16], require that the drift and diffusion coefficients satisfy the global Lipschitz and linear growth assumptions, which are not satisfied by the Heston model with stochastic interest rates. Research on developing time-discrete schemes for the Heston model with stochastic interest rates is scarce. Cozma, Mariapragassam and Reisinger [6] proposed a different log-Euler scheme for the stochastic-local volatility model with stochastic rates and they demonstrated strong convergence without providing a rate. This is the only reference we are aware of concerning the convergence of Monte Carlo algorithms for the Heston model with stochastic interest rates.

In this article, we develop a semi-exact log-Euler scheme for the Heston model with stochastic interest rates, where the driven Brownian motion for interest rate models is independent of the driven Brownian motions for the Heston component. The scheme is an extension of those in Mickel and Neuenkirch [17] and Zheng [22] [23] for the classical Heston model. Under mild assumptions on interest rate models, we show that the underlying scheme converges with order one in the $L^{2}$ norm. The types of options we consider include those with bounded and Lipschitz continuous payoffs and the digital option with a discontinuous payoff. There are two advantages of the scheme we develop. One is that the convergence rate is free of Feller’s index (except for the digital option), i.e., the convergence rate is valid for the full range of parameters in the Heston component of the Heston model with stochastic interest rates. The other advantage is that the convergence rate is higher than the usual convergence rate (one-half in the $L^{2}$ norm) of the standard Euler scheme under standard assumptions. When a numerical scheme has a convergence rate higher than one-half, it is convenient to combine it into Rhee and Glynn’s unbiased estimators. Our result applies to a large class of models, including the Heston-Hull-While model, the Heston-CIR model and the Heston-Black-Karasinski model, among which the Heston-Hull-White model and the Heston-CIR model are particularly attractive in practical applications. We refer readers to Grzelak and Oosterlee [12] for further discussion.

The remainder of this article is organized as follows. In section 2, we review the Heston model with stochastic interest rates and develop a log-Euler scheme for it. Section 3 reviews the unbiased estimators from Rhee and Glynn [20]. In section 4, we derive the relevant convergence rate under several mild assumptions. Section 5 illustrates numerical results to support our theoretical analysis.

Let $(\Omega,\mathcal{F},(\mathcal{F}_{t})_{t\geq 0},\mathbb{P})$ be a filtered probability space satisfying the usual assumptions. The Heston model with stochastic interest rates is of the form

where $(W^{1}_{t})_{t\geq 0}$ and $(W^{2}_{t})_{t\geq 0}$ are two independent $\mathcal{F}_{t}$-adapted Brownian motions and the parameters $k,\theta,\sigma>0$ and $\rho\in[-1,1]$. The initial values $S_{0},V_{0}>0$. Here, $(r_{t})_{t\geq 0}$ is a stochastic interest rate. If we replace $r_{t}$ in the model with a constant, then it becomes the classical Heston model. The classical interest rate models and their generalizations can be found in Brigo and Mercurio [4]. Among them, a large class of interest rate models can be written as

where $\mu,\phi:[0,T]\times\mathbb{R}\rightarrow\mathbb{R}$ are continuous functions and $(W_{t}^{3})_{t\geq 0}$ is a $\mathcal{F}_{t}$-adapted Brownian motion. We assume that $(W_{t}^{3})_{t\geq 0}$ is independent of $(W_{t}^{1})_{t\geq 0}$ and $(W_{t}^{2})_{t\geq 0}$. Furthermore, we assume that there is a unique solution to the equation of $(r_{t})_{t\geq 0}$ above.

Let $X_{t}=\ln(S_{t})$. By using It$\hat{o}$’s lemma, we have

Then substituting the equation of $V_{t}$ into the equation above, we obtain

Since $(V_{t})_{t\geq 0}$ is independent of $(W^{2}_{t})_{t\geq 0}$, the stochastic integral $\int_{0}^{T}\sqrt{V_{t}}\mathrm{d}W_{t}^{2}$ is normally distributed with mean $0$ and variance $\int_{0}^{T}V_{t}\mathrm{d}t$. Therefore, the solution at any finite time horizon $T>0$ can be written as

where $N$ is a standard normal random variable, that is independent of $(V_{t})_{t\in[0,T]}$. Note that $N$ is also independent of $(r_{t})_{t\in[0,T]}$, because the driving Brownian motion $(W^{3}_{t})_{t\in[0,T]}$ of $(r_{t})_{t\in[0,T]}$ is independent of $(W^{2}_{t})_{t\in[0,T]}$ and $(V^{2}_{t})_{t\in[0,T]}$. There are several integrals in equation (1) to be approximated.

It is known that $V_{t}$ follows a scaled noncentral chi-squared distribution given $V_{u}$ for any $u\in[0,t)$, i.e.,

where $\chi_{d}^{2}(\lambda)$ denotes a noncentral chi-squared random variable with degrees of freedom $d=\frac{4k\theta}{\sigma^{2}}>0$ and noncentrality parameter $\lambda>0$ (see Glasserman [11]). Hence, $V_{t}$ can be sampled exactly. Let $(\hat{r}_{ih})_{i=1,..,T/h}$ be an approximate path of $(r_{t})_{t\in[0,T]}$. It is convenient to approximate

by using the Euler scheme based on step size $h$. We denote by $\hat{X}_{T}^{h}$ the approximated solution of $X_{T}$. Let $\hat{S}_{T}^{h}:=e^{\hat{X}_{T}^{h}}$, and then $\hat{S}_{T}^{h}$ is an approximation of $S_{T}$.

In this section, we review the unbiased estimators introduced in Rhee and Glynn [20]. The prices of many options can be expressed as

where $P$ is the payoff function and $Y\in L^{2}$ (i.e., $\mathbb{E}(Y^{2})<\infty$). To estimate this expectation, Rhee and Glynn [20] proposed an estimator

where $\Delta_{n}=Y_{n}-Y_{n-1}$ and $Y_{n}$, $n\in\mathbb{N}$, approximates $Y$ with step size $T/2^{n}$ and $Y_{-1}=0$. In this article, we let

Here, $\mathcal{N}$ is a nonnegative integer-valued random variable that is independent of $Y_{n}$. This estimator is usually referred to as the coupled sum estimator. Other unbiased estimators can be constructed in a similar way (see Rhee and Glynn [20]).

Suppose that $Y_{n}$ converges to $Y$ in the $L^{2}$ norm as $n\rightarrow\infty$. Theorem 1 of Rhee and Glynn [20] showed that if

then $Z$ is an unbiased estimator of $\mathbb{E}(Y)$ (i.e., $\mathbb{E}(Z)=\mathbb{E}(Y)$) with a finite variance. Furthermore, the average computational time of $Z$ is proportional to

Therefore, if $\mathbb{E}\left[(Y_{n}-Y)^{2}\right]=O(2^{-2np})$ with $p>1/2$ (here, $p$ is the convergence rate in the $L^{2}$ norm), we can easily construct a distribution for $\mathcal{N}$ such that $\mathbb{P}(\mathcal{N}\geq n)=O(2^{-n(p+1/2)})$ to ensure that (2) and (3) are finite. The optimal distribution of $\mathcal{N}$ can be calculated by minimizing the product of the variance and the average computational time of $Z$; see Rhee and Glynn [20] and Zheng, Pan and Wang [24] for further discussion. Hence, it is important to investigate the convergence rate of $\mathbb{E}\left[(Y_{n}-Y)^{2}\right]$.

Recall that the interest rate $(r_{t})_{t\in[0,T]}$ follows the stochastic differential equation

where $(W_{t}^{3})_{t\in[0,T]}$ is independent of $(W_{t}^{1})_{t\in[0,T]}$ and $(W_{t}^{2})_{t\in[0,T]}$. Let $c$ and $c_{q}$ be constants regardless of their values, where $c_{q}$ relies on $q\geq 1$. Our analysis throughout this article is based on the following assumption:

Then we proceed to the convergence rate of the underlying log-Euler scheme to approximate the price of an option. Different types of options may have different payoff functions, which can be continuous or discontinuous.