2.1. A new stochastic interest rate model

We now introduce the newly proposed stochastic interest rate model dynamics, which is formulated by making the long-run interest rate level of the Vasicek model another random variable following another OU process. Before that, the Vasicek interest rate model will be presented first, since our model can be viewed as one of its modified versions, and its risk-neutral dynamics under measure Q are given by

$$ \begin{align*} \frac{dS_t}{S_t}&=r_t\,dt+\sigma \,dW_{t}^{S}, \\ dr_t&=\alpha(\beta-r_t)\,dt+\eta \,dW_{t}^{r}, \end{align*} $$

where $S_t$ and $r_t$ are used to denote the underlying price and interest rate, respectively; $\sigma $ and $\eta $ are the constant volatility of the underlying and that of the interest rate, respectively; while $\alpha $ and $\beta $ are the so-called speed of mean-reversion and long-run interest rate level, respectively. Note that $W_{t}^{S}$ and $W_{t}^{r}$ , being two Wiener processes [Reference Wilmott48], are correlated with each other such that $dW_{t}^{S}dW_{t}^{r}=\rho \, dt$ with $\rho $ being the correlation coefficient. As mentioned earlier, we propose to incorporate a stochastic long-run interest level into the Vasicek model so that $\beta $ becomes $\beta _t$ instead of being a constant, and the dynamics of the new model are

$$ \begin{align*} \frac{dS_t}{S_t}&=r_t\,dt+\sigma \,dW_{t}^{S}, \\ dr_t&=\alpha(\beta_t-r_t)\,dt+\eta \,dW_{t}^{r},\\ d\beta_t&=k(\theta-\beta_t)\,dt+\lambda \,dW_{t}^{\beta}. \end{align*} $$

Clearly, the long-run interest level is modelled by an OU process, with its mean-reversion speed, long-run level and volatility being respectively represented by k, $\theta $ and $\lambda $ . We assume that the Wiener process $W_{t}^{\beta }$ does not depend on the other two. To facilitate the computation in the next subsection, the model dynamics that we have just established can be reformulated in a more systematic way as

(2.1) $$ \begin{align} \displaystyle \left[ {\begin{array}{*{20}c} \displaystyle \frac{dS_t}{S_t} \\[.2cm] \displaystyle dr_t \\ \displaystyle d\beta_t\\ \end{array}} \right] = \mu^{Q}\, dt+\Sigma\times C\times\left[ {\begin{array}{*{20}c} dW_{1,t} \\ dW_{2,t} \\ dW_{3,t}\\ \end{array}} \right]. \end{align} $$

Here, we introduce three independent Wiener processes, $W_{1,t}$ , $W_{2,t}$ and $W_{3,t}$ . Here, $\mu ^{Q}$ and $\Sigma $ have the form

$$ \begin{align*} \displaystyle \mu^{Q}=\left[ {\begin{array}{*{20}c} r_t \\ \alpha(\beta-r_t) \\ k(\theta-\beta_t)\\ \end{array}} \right], ~~\Sigma=\left[ {\begin{array}{*{20}c} \sigma & 0 & 0 \\ 0 & \eta & 0 \\ 0 & 0 & \lambda\\ \end{array}} \right], \nonumber \end{align*} $$

and C can be expressed as

$$ \begin{align*} \displaystyle C=\left[ {\begin{array}{*{20}c} 1 & 0 & 0 \\ \rho & \sqrt{1-\rho^2} & 0 \\ 0 & 0 & 1\\ \end{array}} \right] , \end{align*} $$

with which the correlation matrix can be calculated as

$$ \begin{align*}\displaystyle CC^T=\left[ {\begin{array}{*{20}c} 1 & \rho & 0 \\ \rho & 1 & 0 \\ 0 & 0 & 1\\ \end{array}} \right]. \end{align*} $$

After establishing the new stochastic interest rate model, we are going to show how European options can be analytically evaluated under this particular model, the details of which are presented below.

2.2. A closed-form solution

We begin by expressing $U(S,r,\beta ,t)$ , the European option price, by taking the expected value of the option return discounted to the current time

(2.2) $$ \begin{align} U(S,r,\beta,t)=E^{Q}[e^{-\int_{t}^{T}r(s)\,ds}\max(S_T-K,0)|S_t], \end{align} $$

where K denotes the delivery price of the option. The pricing formula presented in (2.2) can be alternatively expressed as

(2.3) $$ \begin{align} U(S,r,\beta,t)=P(r,\beta,t,T)E^{Q^{T}}[\max(S_T-K,0)|S_t], \end{align} $$

if the T-forward measure, $\mathbb {Q}^T$ , is introduced [Reference Brigo and Mercurio4], with $P(r,\beta ,t,T)$ being the value of the bond expiring at time T paying no coupons under the risk-neutral measure. Clearly, in this case, apart from the expectation involved in (2.3), $P(r,\beta ,t,T)$ is also an unknown function that needs to be determined before the option pricing formula can be finally obtained, and its expression is presented in the following theorem.

Theorem 2.1. With the stochastic interest rate following the two-factor model as specified in (2.1), we can compute the value of the bond expiring at time T paying no coupons through

(2.4) $$ \begin{align} \displaystyle P(r,\beta,t,T)=e^{C(\tau)+D(\tau)r_t+E(\tau)\beta_t}, \end{align} $$

where

$$ \begin{align*} \displaystyle C(\tau)&=p_1(\tau)+p_2(\tau)+p_3(\tau)+p_4(\tau)+p_5(\tau)+p_6(\tau),\\ D(\tau)&=-\frac{1}{\alpha}(1-e^{-\alpha\tau}),\\ E(\tau)&=-\frac{\alpha}{k(k-\alpha)}e^{-k\tau}+\frac{1}{k-\alpha}e^{-\alpha\tau}-\frac{1}{k},\\ p_1(\tau)&=\frac{1}{4\alpha}\bigg[\frac{\eta^2}{\alpha^2}+\frac{\lambda^2}{(k-\alpha)^2}\bigg](1-e^{-2\alpha\tau}),\quad p_2(\tau)=\frac{\alpha^2\lambda^2}{4k^3(k-\alpha)^2}(1-e^{-2k\tau}),\\ p_3(\tau)&=-\frac{\alpha\lambda^2}{k(k-\alpha)^2(k+\alpha)}[1-e^{-(k+\alpha)\tau}], \\ p_4(\tau)&=\frac{1}{\alpha}\bigg[-\frac{\eta^2}{\alpha^2}-\frac{\lambda^2}{k(k-\alpha)}+\frac{k\theta}{k-\alpha}\bigg](1-e^{-\alpha\tau}),\\ p_5(\tau)&=\frac{1}{k}\bigg[\frac{\alpha\lambda^2}{k^2(k-\alpha)}-\frac{\alpha\theta}{k-\alpha}\bigg](1-e^{-k\tau}),\quad p_6(\tau)=\bigg(\frac{\eta^2}{2\alpha^2}+\frac{\lambda^2}{2k^2}-\theta\bigg)\tau, \end{align*} $$

with $\tau =T-t$ .

Proof of Theorem 2.1 is provided in Appendix A.

With $P(r,\beta ,t,T)$ being successfully worked out, the task left is to evaluate the expectation in (2.3). Note that the target expectation is considered under $\mathbb {Q}^T$ , the T-forward measure, and thus the model under this particular measure needs to be figured out first, by making use of the dynamics (2.1) under the measure $\mathbb {Q}$ . To conduct measure transformation, it is necessary to find out the numéraires under the two measures. On one hand, $N_{1,t}=e^{\int _{0}^{t}r(s)\,ds}$ is actually the numéraire used under $\mathbb {Q}$ , yielding

$$\begin{align*}\displaystyle \displaystyle dN_{1,t}=N_{1,t}r(t)\,dt,\end{align*}$$

from which one can easily deduce that $\sigma ^{N_{1,t}}=(0,0,0)^T$ is the volatility of $N_{1,t}$ . On the other hand, the numéraire under $\mathbb {Q}^T$ is $\displaystyle N_{2,t}=P(r,\beta ,t,T)$ , which further gives

$$ \begin{align*} dN_{2,t}&=N_{2,t}\bigg[\bigg\{\frac{dC}{dt}+ \frac{dD}{dt}r+\frac{dE}{dt}\beta+\alpha(\beta-r)D+k(\theta-\beta)E+\frac{1}{2}\eta^2D^2+\frac{1}{2}\lambda^2E^2\bigg\}\,dt \\ &\quad +\eta D\,dW_t^r+\lambda E\,dW_t^{\beta}\bigg]. \end{align*} $$

As a result, the volatility term of $N_{2,t}$ is $\sigma ^{N_{2,t}}=(0,\eta D,\lambda E)^T$ . Applying the techniques illustrated in [Reference Brigo and Mercurio4] would certainly lead to the drift term under $\mathbb {Q}^T$ being expressed as

$$ \begin{align*} \mu^{Q^T}&=\mu^{Q}-\Sigma\times CC^T\times\bigg(\frac{\sigma^{N_{1,t}}}{N_{1,t}}-\frac{\sigma^{N_{2,t}}}{N_{2,t}}\bigg) \\ &=\left[ {\begin{array}{*{20}c} r_t \\ \alpha(\beta-r_t)\\ k(\theta-\beta_t) \\ \end{array}} \right]-\left[ {\begin{array}{*{20}c} \sigma & 0 & 0 \\ 0 & \eta & 0 \\ 0 & 0 & \lambda\\ \end{array}} \right]\times\left[ {\begin{array}{*{20}c} 1 & \rho & 0 \\ \rho & 1 & 0 \\ 0 & 0 & 1\\ \end{array}} \right]\times\left[ {\begin{array}{*{20}c} 0 \\ -\eta D(\tau) \\ -\lambda E(\tau)\\ \end{array}} \right] \\ &=\left[ {\begin{array}{*{20}c} r_t-\rho\sigma\eta D(\tau) \\ \alpha(\beta-r_t)+\eta^2D(\tau)\\ k(\theta-\beta_t)+\lambda^2E(\tau)\\ \end{array}} \right]. \end{align*} $$

In this case, we are able to express our model dynamics under $\mathbb {Q}^T$ through

(2.5) $$ \begin{align} \displaystyle \left[ {\begin{array}{*{20}c} \displaystyle \frac{dS_t}{S_t} \\ \displaystyle dv_t \\ \displaystyle dr_t\\ \end{array}} \right] =\left[ {\begin{array}{*{20}c} r_t-\rho\sigma\eta D(\tau) \\ \alpha(\beta-r_t)+\eta^2D(\tau) \\ k(\theta-\beta_t)+\lambda^2E(\tau)\\ \end{array}} \right]dt +\Sigma\times C\times\left[ {\begin{array}{*{20}c} dW_{1,t}^{Q^T} \\[.2cm] dW_{2,t}^{Q^T} \\[.2cm] dW_{3,t}^{Q^T} \\ \end{array}} \right]. \end{align} $$

With the model dynamics being presented in (2.5), the option pricing formula (2.3) can be rearranged as

(2.6) $$ \begin{align} \displaystyle U(y,r,\beta,t)=P(r,\beta,t,T)[P_1-KP_2]. \end{align} $$

Here,

$$ \begin{align*}P_1=\int_{\ln K}^{+\infty}e^{y_T}p(y_T)\,dy_T \quad \text{and}\quad P_2=\int_{\ln K}^{+\infty}p(y_T)\, dy_T\end{align*} $$

with $\displaystyle y_t=\ln (S_t)$ and $p(y_T)$ denoting the $y_T$ density when the current time is t. By denoting the characteristic function as $f(\phi ;t,T,y,r,\beta )$ , we find that

$$ \begin{align*}\displaystyle \frac{e^{y_T}p(y_T)}{f(-j;t,T,y,r,\beta)} \quad\text{with } j=\sqrt{-1}\end{align*} $$

is also a density function, which results in the identity

$$\begin{align*}\displaystyle \int_{-\infty}^{+\infty}e^{y_T}p(y_T)\,dy_T=f(-j;t,T,y,r,\beta).\end{align*}$$

Therefore, we can directly formulate $P_1$ using the connection between the characteristic function and the density through

(2.7) $$ \begin{align} \displaystyle P_1=f(-j;t,T,y,r,\beta)\bigg\{\frac{1}{2}+\frac{1}{\pi}\int_{0}^{+\infty} \text{Real}\bigg[\frac{e^{-j\phi\ln K}f(\phi-j;t,T,y,r,\beta)}{j\phi f(-j;t,T,y,r,\beta)}\bigg]\,d\phi\bigg\}. \end{align} $$

Similarly,

(2.8) $$ \begin{align} \displaystyle P_2=\frac{1}{2}+\frac{1}{\pi}\int_{0}^{+\infty}\text{Real}\bigg[\frac{e^{-j\phi\ln K}f(\phi;t,T,y,r,\beta)}{j\phi}\bigg]\,d\phi. \end{align} $$

Obviously, European call prices under the newly proposed model can be calculated through the derived equations (2.6)–(2.8), whose solution is provided below.

Theorem 2.2. With (2.5) providing the dynamics of the underlying price process, we can compute the objective characteristic function with

(2.9) $$ \begin{align} \displaystyle f=e^{\bar{C}(\phi;\tau)+\bar{D}(\phi;\tau)r_t+\bar{E}(\phi;\tau)\beta_t+j\phi y_t}, \end{align} $$

where

$$ \begin{align*} \displaystyle \bar{D}&=-j\phi D(\tau),\\ \bar{E}&=-j\phi E(\tau),\\ \bar{C}&=-(\phi^2+2j\phi) [p_1(\tau)+p_2(\tau)+p_3(\tau)]+p_7(\tau)+p_8(\tau)+p_9(\tau),\\ p_7(\tau)&=\frac{1}{\alpha}\bigg[(\phi^2+2j\phi)\bigg\{\frac{\eta^2}{\alpha^2}+\frac{\lambda^2}{k(k-\alpha)}\bigg\}+\frac{(\phi^2+j\phi) \rho\sigma\eta}{\alpha}-\frac{k\theta j\phi}{k-\alpha}\bigg](1-e^{-\alpha\tau}),\\ p_8(\tau)&=\frac{1}{k}\bigg[-(\phi^2+2j\phi)\frac{\alpha\lambda^2}{k^2(k-\alpha)}+\frac{\alpha\theta j\phi}{k-\alpha}\bigg](1-e^{-k\tau}),\\ p_9(\tau)&=\bigg[-(\phi^2+2j\phi)\bigg(\frac{\eta^2}{2\alpha^2}+\frac{\lambda^2}{2k^2}\bigg)-\frac{(\phi^2+j\phi)\rho\sigma\eta}{\alpha}+\theta j\phi-\frac{1}{2}\sigma^2(\phi^2+j\phi)\bigg]\tau. \end{align*} $$

We prove Theorem 2.2 in Appendix B.

Till now, one can easily find that European options under the newly constructed stochastic interest rate model can be evaluated with the presented closed-form exact solution. Before the new formula is applied in practice, it is necessary to address two important issues in advance: the first is whether the formula is correct, and the second is how the newly introduced stochastic long-run interest level influences option prices. These results are illustrated in the next section.

Figure 1 Our price versus Monte Carlo price.

4.1. Data description

S&P 500 index options within the period of January–June 2018 are selected for carrying out empirical analysis. For parameter estimation purposes, it is necessary to first figure out what is the price of an option, which is actually selected as the mean of ask and bid prices, by conversion. The U.S. Treasury Bill Rate of three months that is announced every day is selected as the value of the current interest rate [Reference Shu and Zhang43]. Moreover, before the parameters are estimated, we need to apply several filters on the raw data so that market data noise can be eliminated [Reference Bakshi, Cao and Chen2].

First, only one-day data of a week are used for parameter estimation, because calibrating option pricing models can be time costly and this would enable a longer period to be studied. In particular, we choose Wednesday options data since among the five working days of each week, there is a least probability that Wednesday is a holiday and the day-of-the-week effect can be possibly eliminated. Thursday options data are also adopted, but for the purpose of assessing the model performance in the out-of-sample manner, that is, comparing the Thursday prices calculated using the parameters determined for Wednesday and those listed in real markets. Second, options far from the maturity (more than 90 days) and with less than 7 days to maturity are both discarded because the former category is not popular in real markets due to their high prices while the latter has less time value. Third, as a result of the illiquidity associated with very deep out-of-money and in-the-money options, options with more than 10% moneynessFootnote ² and less than $-$ 10% moneyness are all deleted. At last, due to the price instability, we remove options whose values are less than $1/8.

Having applied all the filters to the raw data set, we can now proceed to determine model parameters through the filtered data, and the related issues are discussed below.

4.2. Parameter estimation

Parameter estimation is always the first and also one of the most important steps in assessing model performance. One common approach in estimating model parameters is to determine a set of “optimal” parameters, with which the produced model prices are the “closest” to corresponding market option prices. Thus, to measure the “closeness”, we actually need an appropriate function to calculate the “distance” between market and model prices. Following [Reference Christoffersen and Jacobs10], dollar mean-squared errors are selected for such a distance, the expression of which can be written as

(4.1) $$ \begin{align} \displaystyle MSE=\frac{1}{N}\sum_{i=1}^{N}(C_{\text{Market}}-C_{\text{Model}})^2, \end{align} $$

where $C_{\text {Model}}$ and $C_{\text {Market}}$ respectively represent the option price produced by the model and that listed on the market, and N denotes the number of options used when estimating the parameters of that day.

Clearly, the parameter determination problem reduces to finding the optimal parameters that can lead to the objective function (4.1) to be minimized, solving which requires the selection of an appropriate optimization approach. In particular, local minimization algorithms are not appropriate, since the convexity of (4.1) cannot be guaranteed and the existence of different local minima is possible. To avoid ending up with a local minimum, global optimization algorithms are much more preferred.

What we choose is adaptive simulated annealing (ASA) [Reference Ingber31], an improved version of simulated annealing. The history of this algorithm dates back to 1983, when simulated annealing was initially established for highly nonlinear problems [Reference Kirkpatrick, Gelatt and Vecchi33]. It was then improved in [Reference Szu and Hartley45] by developing fast simulated annealing, an algorithm that ensures attainability of the global minimum in a finite time period. A further improvement was made in [Reference Ingber28] by proposing ASA, formerly known as “very fast simulated reannealing”, and it is able to adjust random step selection automatically based on the actual progress of the algorithm so that the speed of the algorithm can further be accelerated. In fact, ASA has already been implemented in different areas [Reference Chen, Luk and Liu8, Reference Chen and Luk9], including the application in calibrating models used to price options [Reference Ingber29, Reference Mikhailov and Nögel40]. Furthermore, the open-source code listed on Ingber’s homepage [Reference Ingber31] is updated regularly according to user feedback so that its flexibility and powerfulness are promised.

Moreover, as mentioned earlier, we also want to emphasize the mean reversion property of the stochastic long-term interest rate. Thus, we construct another two-factor Vasicek model using a normal distribution, and combine it with the Black–Scholes model (the BSVN model hereafter) as another benchmark:

$$ \begin{align*} \frac{dS_t}{S_t}&=r_t\,dt+\sigma \,dW_{t}^{S}, \\ dr_t &=\alpha(\beta_t-r_t)\,dt+\eta \,dW_{t}^{r},\\ d\beta_t&=\bar{k}\,dt+\lambda \,dW_{t}^{\beta}. \end{align*} $$

By implementing ASA, we report the extracted parameters (averaged daily) in Table 1 Footnote ³.

Table 1 Estimated parameters.

4.3. Empirical comparison

The assessment of the model performance based on the parameters estimated from real market data is now provided. If the difference between market and model prices is lower, we can regard the performance of that model better. There are mainly two categories in defining such pricing differences, yielding in-sample and out-of-sample errors. In particular, the difference between model prices produced using parameters estimated and market prices used for parameter estimation are defined as in-sample errors, to serve a guide of whether the target model can produce a good market fitness. In contrast, the difference between option prices computed with parameters determined from other sets of data and market option prices that have not been used to estimate parameters are defined as out-of-sample errors, to check whether the model can provide consistent results across different data sets so as to assess its prediction ability.

We compare in-sample as well as out-of-sample errors of the BSV model, the BSVN model and our model in Table 2. Note that the performance of the BSV model is generally worse than that of our model, in the sense that both errors are lower under our model, with a significant 20% improvement in the in-sample errors. Moreover, although the in-sample fitness of the BSVN model has been improved by the introduction of an additional normal factor, the out-of-sample performance is much worse than the other two models. This is actually strong evidence supporting the use of a mean reverting process to model the stochastic long-term interest rate. It should also be pointed out that the magnitude of out-of-sample errors are higher than that of in-sample errors under both models, which is reasonable since the optimal parameters obtained from one data set are not necessarily optimal for another set.

Table 2 In- and out-of-sample errors.

We perform a t-test to see if the introduction of a two-factor stochastic interest rate and the use of a mean-reversion process in modelling the long-run equilibrium level can significantly reduce the pricing errors. The corresponding results are provided in Table 3. One can clearly see that the in- and out-of-sample errors differ significantly between our model and the BSV model at 1% and 4% significance levels, respectively, while those between our model and the BSVN model also differ significantly at 7% and 3% significance levels, respectively. This demonstrates the overall better performance of our model over both the BSV and BSVN models.

Table 3 Significance tests of pricing errors.

Table 4 Out-of-sample errors across different moneyness.

Figure 4 Implied volatility.

One may also be interested in how the three models perform differently when out-of-sample errors are classified into different moneyness. This is because a series of strike prices are available even for traded options of a day. We abbreviate “out of money”, “at the money” and “in the money” in the parentheses by (O), (A) and (I), respectively, and the results are presented in Table 4. Although the out-of-sample performance of the BSV model is slightly better than that of our model, as far as out of money options are concerned, the BSV model provides much worse performance in the other two categories compared with our model, with the greatest improvement shown by the at the money options. However, the out-of-sample performance of the BSVN model is the worst across all the three categories, which again demonstrate the importance of mean reversion in the stochastic long-term interest rate.

Considering that the accurate determination of implied volatility is a vital issue in finance practice, we display the ability of the three models in reproducing market implied volatility in Figure 4. A clear pattern in our model provides a better match to market-implied volatility, which emphasizes the significance of considering stochastic long-term interest rate as well as its mean reversion property.