A Stochastically Correlated Bivariate Square-root Model

1. Introduction

Scalar square-root diffusion processes, also called CIR processes, due to [1], have two properties that are important in applications in finance: (i) they ensure mean reversion of the state variables towards a long run level and (ii) unlike the [2] model, they do not take negative values. The importance of such properties clearly replicates in the multi-dimensional scenario, especially when correlation is allowed among the entries of the process.

In turn, the stochastic models found in finance, especially in interest rate markets, usually assume uncorrelated state variables contradicting both intuition and historical data. In fact, considerable difficulties arise in the pricing procedures when entering with correlation in the model, even when numerical solutions of the corresponding partial differential equations are given by finite difference or finite element methods. Concerning this matter, we quote [3]: “requiring $d{W}_{1}d{W}_2=\rho dt$ with $\rho \ne 0$ in the two-factor CIR model would indeed destroy analytical tractability: It would no longer be possible to compute analytically bond prices and rates starting from the short rate factors.” On the other hand, as it is still argued by [3], whenever the correlation plays a more relevant role or when higher precision is needed in the pricing procedure, we need to resort to models that allow for more realistic patterns.

With this in mind, we present a bi-dimensional (two-factor) mean-reversion stochastically correlated process in which entries are produced by two concatenated CIR-type models. Such process is affine according to [4].

We obtained the closed form solution - of exponential format - of the conditional moment generating function and the conditional characteristic function of the density function of the integral of the process governed by the model as well as those of the process itself at final time T. A key point that permitted us to apply the above transforms in the problem of pricing derivatives, particularly with respect to fixed income markets, was devising that the integral of the process - and not the process per se - was the adequate mathematical object to achieve the results.

It is noteworthy that by setting our stochastic correlation to zero we reduce our model to a pair of scalar uncorrelated square-root (CIR) diffusions. We believe that the proposed model encompasses a range of processes similar to that produced by a bivariate CIR process in which entries are correlated in the standard way, i.e., exhibiting a constant correlation coefficient, with the advantage of exhibiting an analytical solution.

We applied our results in the context of the interest rate financial markets where the interest rate process stems from the two-factor stochastically correlated CIR model we present in the paper. We obtained closed-form expressions for the price of a zero-coupon bond and, in the sequel, for the price of the IDI option, an Asian option that commonly appears in the Brazilian fixed-income markets. For the latter case, we used the fast and accurate pricing method based on Fourier series expansions proposed by Fang and Osterlee [5] and applied it to interest rate options in [6].

Via numerical simulations, we showed the (good) behavior of the two-factor process in the face of a certain stochastic correlation condition, under stringent circumstances. Moreover, a result concerning the effect of the sensitivity parameter ($ϵ$) on the stochastic correlation coefficient ($\varrho$), and ultimately, the possibility of emulating bivariate CIR processes with constant correlation coefficient, was unveiled. Focusing on the financial scenario, we analyze the impact of the correlation parameter ($\rho$) on the conditional probability density functions associated with the interest rate process - which in turn stems from the two-factor process - obtaining again a good confirmation (non-negativity). Finally, some option price curves are obtained, parameterized by $\rho$. We also used real data from the financial market in the model design.

The mathematical tools we relied on are the Feynman-Kac formula, the separation of variables approach, the Riccati equations with known explicit solutions and the series expansions proposed by [5].

The rest of this paper is organized as follows: Section 2 presents the stochastically correlated two-factor affine model and its probability background. In Section 3 we derive the expressions of the conditional moment generating function and the conditional characteristic function associated to the model. Section 4 addresses applications to finance, and Section 5 the numerical results.

2. The Model

Let $(\Omega ,\mathcal{F},\mathbb{Q},{\left\{{\mathcal{F}}_t\right\}}_{t\in [0,T]})$ be a filtered space where ${\left\{{\mathcal{F}}_t\right\}}_{t\in [0,T]}$ is a filtration, $\mathbb{Q}$ is a probability measure, and the adapted processes $Z_1=(Z_1\left(t\right),t\in [0,T])$ and $Z_2=(Z_2\left(t\right),t\in [0,T])$ are independent standard Brownian Motions under $\mathbb{Q}$ (see, e.g., [7]).

We introduce the pair of square-root diffusions, in which the underlying variables are stochastically correlated, as follows. For $t\in [0,T]$ and initial conditions $x_1\left(0\right)={\overline{x}}_1$ and $x_2\left(0\right)={\overline{x}}_2$,

$$\begin{array}{ccc}\hfill d{x}_1\left(t\right)& =& {\kappa }_1({\theta }_1-x_1\left(t\right))dt+{\sigma }_1\sqrt{x_1\left(t\right)}d{W}_1\left(t\right)\hfill \end{array}$$

(1)

and

$$\begin{array}{ccc}\hfill d{x}_2\left(t\right)& =& {\kappa }_2({\theta }_2-x_2\left(t\right))dt+{\sigma }_2\sqrt{(1-{\rho }^2)}\sqrt{x_2\left(t\right)}d{W}_2\left(t\right),\hfill \end{array}$$

(2)

where ${\overline{x}}_1$ and ${\overline{x}}_2$ are prescribed random variables and $\rho \in [0,1)$. We have that ${\kappa }_i$, ${\theta }_i$ and ${\sigma }_i$ are positive real numbers. Each stochastic process $x_i\left(t\right)$ is pulled to a level ${\theta }_i$ at a rate ${\kappa }_i$ with a volatility ${\sigma }_i$. We have that $W_1$ and $W_2$ are stochastically correlated processes, where $W_1\left(t\right)=Z_1\left(t\right)$ is a Brownian Motion and

$$\begin{array}{c}\hfill W_2\left(t\right)={\int }_0^t\varrho \left(t\right)d{Z}_1\left(t\right)+{\int }_0^t\sqrt{1-\varrho {\left(t\right)}^2}d{Z}_2\left(t\right)\end{array}$$

(3)

is a Itô process.

Hence,

$$\begin{array}{c}\hfill d[W_1,W_1]\left(t\right)=d[W_2,W_2]\left(t\right)=dt,\mathrm{and}d[W_1,W_2]\left(t\right)=\varrho \left(t\right)dt.\end{array}$$

(4)

The correlation $\varrho$ is specified via the adapted stochastic process

$$\begin{array}{c}\hfill \varrho \left(t\right)=\frac{\rho ϵ}{\sqrt{x_1\left(t\right)x_2\left(t\right)(1-{\rho }^2)}},\end{array}$$

(5)

$\rho \in (-1,1)$. Note that if $\rho =0$, then $W_2=Z_2$ and we recover the standard uncorrelated case. If we change the specification of $\varrho$ to a constant, say, $\varrho =\rho$, then $W_2=\rho Z_1+\sqrt{1-{\rho }^2)}{Z}_2$ and $E\left[W_2^2\left(t\right)\right]=t$, so $W_2$ becomes a standard Brownian Motion.

The dynamics (1) & (2) also reads

$$\left[\begin{array}{c}d{x}_1\left(t\right)\\ d{x}_2\left(t\right)\end{array}\right]=\left[\begin{array}{c}{\kappa }_1({\theta }_1-x_1\left(t\right))\\ {\kappa }_2({\theta }_2-x_2\left(t\right))\end{array}\right]dt+\overline{\sigma }(x_1\left(t\right),x_2\left(t\right))\left[\begin{array}{c}d{Z}_1\left(t\right)\\ d{Z}_2\left(t\right)\end{array}\right],$$

(6)

where

$$\overline{\sigma }(x_1\left(t\right),x_2\left(t\right))=\left[\begin{array}{cc}{\sigma }_1\sqrt{x_1\left(t\right)}& 0\\ 0& {\sigma }_2\sqrt{(1-{\rho }^2)}\sqrt{x_2\left(t\right)}\end{array}\right]\left[\begin{array}{cc}1& 0\\ \varrho \left(t\right)& \sqrt{1-{\left(\varrho \left(t\right)\right)}^2}\end{array}\right]$$

(7)

and $Z_1$ and $Z_2$ are independent standard Brownian Motions, as defined at the beginning of the section. Using the specification (5) and after some algebraic manipulation, it follows that $\overline{\sigma }(x_1\left(t\right),x_2\left(t\right))\overline{\sigma }{(x_1\left(t\right),x_2\left(t\right))}^{\prime }$ is of affine form.

In the above equations, ${\kappa }_i$ are the speeds of adjustment, ${\theta }_i$ are the long-term means, ${\sigma }_i$ are the volatility factors, $\rho$ is a correlation parameter and $ϵ$ is a sensitivity factor - a useful parameter as we shall see in Section 5. Equations (1) and (2) produce Feller processes known in the context of finance as the CIR processes due to [1]. According to [8], the following condition must hold in order to ensure a strong solution for (1) and (2) and to avoid the possibility of negative values of $x_1$ and $x_2$.

$$\begin{array}{c}\hfill \frac{2{\kappa }_1{\theta }_1}{{\sigma }_1^2}+\frac{2{\kappa }_2{\theta }_2}{{\sigma }_2^2(1-{\rho }^2)}\ge 1.\end{array}$$

(8)

Moreover, we rely on having

$$\begin{array}{c}\hfill \varrho {\left(t\right)}^2\le 1a.s.\forall t\in [0,T].\end{array}$$

(9)

In order to conform with the condition above, we benefit from the existence of the sensitivity parameter $ϵ$, setting it as follows.

$$ϵ\in \left[0,(1-{\rho }^2)\sqrt{min({\overline{x}}_1,{\theta }_1)}\sqrt{min({\overline{x}}_2,{\theta }_2)}\right].$$

(10)

The comprehensive numerical results of Section 5 assert that the above condition suffices to ensure the good behavior of the correlation process and that the bivariate process x covers a broad range of CIR-type correlated processes.

3. Underpinning Results

The theorems that follow consider the stochastically correlated two-factor process $x\left(t\right)={(x_1\left(t\right),x_2\left(t\right))}^{\prime }$, as given above, relying on the fact that this process is affine. Here ’ stands for the transpose.

The conditional moment generating function and the characteristic function of Theorem 1 will be applied in Section 4 in order to calculate, respectively, the price of a zero coupon bond and the price of an interest rate option, namely, the IDI option. In both cases, the interest rate is given by the additive process (29), stemming from the two-factor process x.

Theorem 1. 

The closed-form solution of the conditional characteristic function (resp. conditional moment generating function) associated with the $[t,T]$-integral of the two-factor process x with coefficient $\eta ={({\eta }_1,{\eta }_2)}^{\prime }\in {\mathbb{R}}^2$, namely ${\int }_t^T{\eta }^{\prime }x\left(s\right)ds$ (resp. $-{\int }_t^T{\eta }^{\prime }x\left(s\right)ds$), is given by

$$\begin{array}{c}\hfill \widehat{f}(x\left(t\right),t,u)=\mathbb{E}[e^{iu{\int }_t^T{\eta }^{\prime }x\left(s\right)ds}\mid x\left(t\right)]=e^{\alpha (t,T)-{\beta }_1(t,T)x_1\left(t\right)-{\beta }_2(t,T)x_2\left(t\right)}\end{array}$$

(11)

$$\begin{array}{c}\hfill (resp.\\ \hfill G(x\left(t\right),t,u)=\mathbb{E}[e^{-u{\int }_t^T{\eta }^{\prime }x\left(s\right)ds}\mid x\left(t\right)]=e^{\alpha (t,T)-{\beta }_1(t,T)x_1\left(t\right)-{\beta }_2(t,T)x_2\left(t\right)})\end{array}$$

(12)

where

$$\begin{array}{ccc}\hfill \alpha (t,T)& =& {\int }_t^T-{\kappa }_1{\theta }_1{\beta }_1(t,T)-{\kappa }_2{\theta }_2{\beta }_2(t,T)+\rho ϵ{\sigma }_1{\sigma }_2{\beta }_1(t,T){\beta }_2(t,T)ds,\hfill \end{array}$$

(13)

$$\begin{array}{ccc}\hfill {\beta }_j(t,T)& =& \frac{2{\eta }_j\left(e^{{\gamma }_j(T-t)}-1\right)}{({\gamma }_j+{\kappa }_j)e^{{\gamma }_j(T-t)}+({\gamma }_j-{\kappa }_j)},j=1,2,\hfill \end{array}$$

(14)

with

$$\begin{array}{ccc}\hfill {\gamma }_1& =& \sqrt{{\kappa }_1^2+2i{\sigma }_{1}u{\eta }_1}and{\gamma }_2=\sqrt{{\kappa }_2^2+2i{\sigma }_2(1-{\rho }^2)u{\eta }_2}\hfill \end{array}$$

(15)

$$\begin{array}{ccc}\hfill (resp.{\gamma }_1& =& \sqrt{{\kappa }_1^2+2{\sigma }_{1}u{\eta }_1}and{\gamma }_2=\sqrt{{\kappa }_2^2+2{\sigma }_2(1-{\rho }^2)u{\eta }_2}).\hfill \end{array}$$

(16)

Proof. 

We address the proof for the moment generating function case. The proof for the characteristic function case follows closely. So, invoking [4], the moment generating function of a multifactor - say m-dimensional, stochastic differential equation (SDE) of the form

$$dx\left(t\right)=\mu \left(x\right(t\left)\right)dt+\sigma \left(x\right(t\left)\right)dW\left(t\right),$$

has an exponential affine solution if the above coefficients are of the form

$$\mu \left(x\right)=K_0+K_{1}x\mathrm{and}{\left(\sigma \left(x\right)\sigma {\left(x\right)}^{\prime }\right)}_{ij}=H_{0ij}+H_{1ij}x,$$

where $K_0$, $K_1$, $H_{0ij}$ and $H_{1ij}$ have adequate dimensions and do not depend on x.

The above SDE captures (1) and (2) setting $m=2$, $K_0={[{\kappa }_1{\theta }_1;{\kappa }_2{\theta }_2]}^{\prime }$, $K_1=\left[\begin{array}{cc}-{\kappa }_1& 0\\ 0& -{\kappa }_2\end{array}\right]$ and

$$\left(\sigma \left(x\right)\sigma {\left(x\right)}^{\prime }\right)=\left[\begin{array}{cc}{\sigma }_1^2{x}_1& \rho {\sigma }_1{\sigma }_2ϵ\\ \rho {\sigma }_1{\sigma }_2ϵ& {\sigma }_2^2\left(1-{\rho }^2\right)x_2\end{array}\right],$$

i.e., the coefficients are of affine form. In order to have $\det \left(\sigma \left(x\right)\sigma {\left(x\right)}^{\prime }\right)>0$ we require having

$$\begin{array}{c}\hfill x_1{x}_2(1-{\rho }^2)>{\rho }^2{ϵ}^2.\end{array}$$

(17)

Now, applying the Feynman-Kac formula to the second expression of (12), leads us to

$$\frac{\partial G\left(x\right(t),t,u)}{\partial t}+\nabla G{(x\left(t\right),t,u)}^{\prime }(K_0+K_{1}x)+\sum _{i,j}\frac{\partial G\left(x\right(t),t,u)}{\partial x_i\partial x_j}(H_{0ij}+H_{1ij}x)=u{\eta }^{\prime }xG(x\left(t\right),t,u),$$

and so, invoking the last expression of Equation (11), we have that

$$\begin{array}{c}\frac{\partial \alpha }{\partial t}+x_1\left(t\right)\frac{\partial {\beta }_1}{\partial t}+x_2\left(t\right)\frac{\partial {\beta }_2}{\partial t}+\left[{\kappa }_1({\theta }_1-x_1\left(t\right))\right]{\beta }_1+\left[{\kappa }_2({\theta }_2-x_2\left(t\right))\right]{\beta }_2+\frac{1}{2}{\sigma }_1^2{x}_1\left(t\right){\beta }_1^2\hfill \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em} +\frac{1}{2}{\sigma }_2^2{x}_2\left(t\right){\beta }_2^2(1-{\rho }^2)+\rho {\sigma }_1{\sigma }_2{\beta }_1{\beta }_2ϵ=u{\eta }^{\prime }x\left(t\right).\hfill \end{array}$$

Gathering the terms in $x_1$ and $x_2$ gives us

$$\frac{\partial \alpha }{\partial t}=-{\kappa }_1{\theta }_1{\beta }_1-{\kappa }_2{\theta }_2{\beta }_2-\rho {\sigma }_1{\sigma }_2{\beta }_1{\beta }_2ϵ$$

(18)

and the Riccati differential equations

$$\frac{\partial {\beta }_1}{\partial t}=-{\kappa }_1{\beta }_1-\frac{{\sigma }_1^2}{2}{\beta }_1^2+{\eta }_1$$

(19)

and

$$\frac{\partial {\beta }_2}{\partial t}=-{\kappa }_2{\beta }_2-\frac{{\sigma }_2^2}{2}{\beta }_2^2(1-{\rho }^2)+{\eta }_2.$$

(20)

Now, for terminal conditions

$$\alpha (T,T)={\beta }_1(T,T)={\beta }_2(T,T)=0$$

After some algebraic manipulation, we obtain that the solutions of the Riccati Equations (19) and (20) are indeed given by Equations (13) and (14). □

The characteristic function version of the following theorem will be applied in Section 5 in order to calculate the probability density function of the additive process given by Equation (29) exhibited ahead.

Theorem 2. 

The closed-form solution of the conditional characteristic function transform associated with the two-factor process x at time T with coefficient $\eta ={({\eta }_1,{\eta }_2)}^{\prime }\in {\mathbb{R}}^2$, namely, ${\eta }^{\prime }{x}_T$, is given by

$$\begin{array}{c}\hfill \widehat{f}(x\left(t\right),t,u)=\mathbb{E}[e^{iu{\eta }^{\prime }x\left(T\right)}\mid x\left(t\right)]=e^{(\alpha (t,T)-{\beta }_1(t,T)x_1\left(t\right)-{\beta }_2(t,T)x_2\left(t\right))},\end{array}$$

(21)

where

$$\begin{array}{ccc}\hfill \alpha (t,T)& =& {\int }_t^T-{\kappa }_1{\theta }_1{\beta }_1\left(s\right)-{\kappa }_2{\theta }_2{\beta }_2\left(s\right)+\rho ϵ{\sigma }_1{\sigma }_2{\beta }_1\left(s\right){\beta }_2\left(s\right)ds,\hfill \end{array}$$

(22)

$$\begin{array}{ccc}\hfill {\beta }_1(t,T)& =& \frac{-2{\kappa }_1{e}^{-{\kappa }_{1}t}iu\eta }{{\sigma }_1^2(e^{-{\kappa }_{1}t}-1)iu\eta -2{\kappa }_1},\hfill \end{array}$$

(23)

$$\begin{array}{ccc}\hfill {\beta }_2(t,T)& =& \frac{-2{\kappa }_2{e}^{-{\kappa }_{2}t}iu\eta }{{\sigma }_2^2(1-{\rho }^2)(e^{-{\kappa }_{2}t}-1)iu\eta -2{\kappa }_2}.\hfill \end{array}$$

(24)

Proof. 

Following the same steps as those of the proof of Theorem 1, we obtain the Riccati differential equations

$$\frac{\partial {\beta }_1}{\partial t}=-{\kappa }_1{\beta }_1-\frac{{\sigma }_1^2}{2}{\beta }_1^2$$

(25)

and

$$\frac{\partial {\beta }_2}{\partial t}=-{\kappa }_2{\beta }_2-\frac{{\sigma }_2^2}{2}{\beta }_2^2(1-{\rho }^2),$$

(26)

and given the terminal conditions

$$\alpha (T,T)=0,\mathrm{and}{\beta }_1(T,T)=iu{\eta }_1,\mathrm{and}{\beta }_2(T,T)=iu{\eta }_2,$$

Again, after some algebraic manipulation, we obtain that Equations (23) and (24), respectively, satisfy the Riccati Equations (25) and (26). □

Remark 1. 

It is worth mentioning that the two-factor uncorrelated CIR process given by

$$\begin{array}{ccc}\hfill d{x}_1\left(t\right)& =& {\kappa }_x({\theta }_1-x_1\left(t\right))dt+{\sigma }_1\sqrt{x_1\left(t\right)}d{W}_1\left(t\right),\hfill \end{array}$$

(27)

$$\begin{array}{ccc}\hfill d{x}_2\left(t\right)& =& {\kappa }_r({\theta }_2-x_2\left(t\right))dt+{\sigma }_2\sqrt{x_2\left(t\right)}d{W}_2\left(t\right),\hfill \end{array}$$

(28)

introduced in the famous paper [1] which closed form solution is shown in [3], is a particular case of the analytical solution (11) with $\rho =0$.

4. Applications to Finance: Fixed Income Markets with Decomposed Nominal Rates

We consider a financial interest rate market underpinned by the short rate stochastic process $R=\left(R\right(t),t\in [0,T\left]\right)$ which stems from the two-factor stochastically correlated (affine) CIR model given by (1)–(4), namely

$$\begin{array}{c}\hfill R\left(t\right)={\eta }_1{x}_1\left(t\right)+{\eta }_2{x}_2\left(t\right),\end{array}$$

(29)

${({\eta }_1,{\eta }_2)}^{\prime }\in {\mathcal{R}}^2$. We assume that $\mathbb{Q}$ is a risk-neutral probability measure.

Now, recalling the Fisher effect [9] in economics, we have that the nominal short rate $R\left(t\right)$ can be decomposed into the real interest rate and the inflation expectation. This conforms to our setting, simply particularizing $({\eta }_1,{\eta }_2)=(1,1)$ and denoting, for convenience, the real interest rate and the inflation by $x_1\left(t\right)$ and $x_2\left(t\right)$ respectively.

Under this framework, and relying on Theorem 1, we shall obtain the price of a bond and an interest rate option as follows.

4.1. Bond Pricing

The zero-coupon bond is an interest rate derivative that pays 1 at the expiration date T and no other payment before T. Relying on the seminal work [10], the no-arbitrage price at time $t\in [0,T]$ of the zero-coupon bond expiring at T is given by

$$\begin{array}{c}\hfill D(t,T)={\mathbb{E}}^{\mathbb{Q}}\left[e^{-{\int }_t^{T}R\left(s\right)ds}|R\left(t\right)\right].\end{array}$$

(30)

Proposition 1. 

Consider an interest rate market with nominal short rate given by Equation (29), $({\eta }_1,{\eta }_2)=(1,1)$, where the real yield $x_1\left(t\right)$ and inflation premium $x_2\left(t\right)$ follows the stochastic correlated square-root diffusion model (1)–(4). Then the dynamic of the no-arbitrage price of the bond at the time interval $[0,T]$ is given by

$$\begin{array}{cc}\hfill & dD(t,T)=D(t,T)\{\left[{\eta }_1{x}_1\left(t\right)+{\eta }_2{x}_2\left(t\right)\right]dt\hfill \\ \hfill & +\left[{\sigma }_1{\beta }_1(t,T)d{W}_1^{\mathbb{Q}}\left(t\right)+{\sigma }_2(1-{\rho }^2){\beta }_2(t,T)d{W}_2^{\mathbb{Q}}\left(t\right)\right]\}.\hfill \end{array}$$

(31)

Proof. 

Using the moment generating function version of Theorem 1, we have that the solution to Equation (30) is of exponential-affine form, as given by the right side of (11) with $({\eta }_1,{\eta }_2)=(1,1)$. Relying on Ito’s formula and after some algebraic manipulation, it follows that

$$\begin{array}{cc}\hfill & dD(t,T)=D(t,T)[-{\alpha }_t(t,T)+{\beta }_1^{\prime }(t,T)x_1\left(t\right)+{\beta }_2^{\prime }(t,T)x_2\left(t\right)\hfill \\ \hfill & -{\kappa }_1({\theta }_1-x_1\left(t\right)){\beta }_1(t,T)-{\kappa }_2({\theta }_2-x_2\left(t\right)){\beta }_2(t,T)+\hfill \\ \hfill & +\frac{1}{2}{\sigma }_1^2{x}_1\left(t\right){\beta }_1{(t,T)}^2+\frac{1}{2}{\sigma }_2^2(1-{\rho }^2)x_2\left(t\right){\beta }_2{(t,T)}^2+{\sigma }_1{\sigma }_2\rho ϵ{\beta }_1(t,T){\beta }_2(t,T)]dt+\hfill \\ \hfill & +D(t,T)\left[{\sigma }_1{\beta }_1(t,T)d{W}_1^{\mathbb{Q}}\left(t\right)+{\sigma }_2(1-{\rho }^2){\beta }_2(t,T)d{W}_2^{\mathbb{Q}}\left(t\right)\right].\hfill \end{array}$$

(32)

Using equation (19) and (20) for ${\beta }_1^{\prime }(t,T)$ and ${\beta }_2^{\prime }(t,T)$, respectively, we obtain Equation (31). □

We may notice that the drift that ultimately appears in dynamic (31) is the short rate process itself, as it ought to be.

As stated by Brigo and Mercurio ([3]), and shown in [11] for example, a two-factor model is capable of explaining up to $91\%$ of the variations in the yield curve. Such a model is needed to provide a realistic evolution of the zero-coupon bond curve. So, we contribute to this, by means of sophisticating the modeling of the term structure of the interest rates via the additive two-factor model given by (1)-(4), (29) and $({\eta }_1,{\eta }_2)=(1,1)$. We also provide a leak concerning the statement made in [3], in which the authors say that “requiring $d{W}_{1}d{W}_2=\rho dt$ with $\rho \ne 0$ in the two-factor CIR model would indeed destroy analytical tractability: It would no longer be possible to compute analytically bond prices and rates starting from the short rate factors.” Indeed, we have shown that bond prices are analytically tractable under the two-factor correlated CIR model mentioned above.

5. Numerical Results

In this section, we discuss some issues of the stochastic correlated square-root diffusion process given by Equations(1)–(4). In Section 5.1, we study the model’s behavior in the face of condition given by the inequality (9) - which involves the stochastic correlation coefficient $\varrho$. In (Section 5.2) we study the variation of the stochastic correlation coefficient $\varrho$ in time and in the probability space for fixed parameters $\rho$ and $\epsilon$. In the sequel we let these parameters vary.

Via Fourier series, we proceed to present the conditional probability density function of the integral of the short rate process R - the weighted sum of the components of the two-factor CIR model (1)–(4) (see Equation (29)). The conditional probability density function is equally provided for R at terminal date T. Such sort of conditional density functions are the cornerstones for obtaining the price of financial derivative contracts. We shall see that, in the long run especially, the correlation parameter $\rho$ plays a sensitive role in the volatility of the zero coupon bond, and therefore on the prices of interest rate derivatives in general. In the sequel, we obtain the price of the IDI option (Section 5.4).

Nonetheless, the important stochastic correlation analysis will be carried out varying $ϵ$ in the whole range given by (33). Our dynamic is as given by Equations (1)–(4) with ${\kappa }_1=0.15$, ${\kappa }_2=0.15$, ${\theta }_1=0.05$, ${\theta }_2=0.05$, ${\sigma }_1=0.015$, ${\sigma }_2=0.025$, $T-t=1$ year, and initial condition $x_1\left(0\right)=0.05$ and $x_2\left(0\right)=0.05$. We also set ${\eta }_1={\eta }_2=1$. We benefit from the existence of the sensitivity parameter $ϵ$, setting it equal to the maximal value of the range given by (10), i.e.,

$$\begin{array}{c}\hfill ϵ=(1-{\rho }^2)\sqrt{min(x_1\left(0\right),{\theta }_1)}\sqrt{min(x_2\left(0\right),{\theta }_2)}.\end{array}$$

(33)

This maximum value will hold throughout the section. Nonetheless, the important stochastic correlation analysis will be carried out varying $ϵ$ in the whole range given by (10).

Figure 1 illustrates a bivariate sample path evolving according to (1)–(4), with $\rho =0.5$ and $ϵ$ as above.

5.1. The Stochastic Correlation Condition (9) is Mild, at Least at a Certain Niche

To verify if the stochastic correlation coefficient $\varrho \left(t\right)$ performs according to condition (9), we produced 100.000 path simulations of the stochastic correlated CIR process with time discretization $\Delta t=0.01$ years and period of one year, using the simple Euler scheme as given in [12].

Firstly, we used [13], which provided us with parameter estimates based on real-world financial data. There, the authors used a three-factor uncorrelated CIR process to model zero-coupon bonds. We tested all three possible combinations of the parameters presented in Table 2 of the referenced paper, not finding a rate over $0.01\%$ concerning pointwise violations of the stochastic correlation condition (9) throughout $[t,T]$, for a variety of values of $\rho$.

In order to test a more restrictive scenario, we set ${\kappa }_1={\kappa }_2=0.75$, ${\theta }_1={\theta }_2=0.01$ and ${\sigma }_1={\sigma }_2=0.03$ for $x_1\left(t\right)={\theta }_1$ and $x_2\left(t\right)={\theta }_2$. While agreeing with condition (8), the setup above enforces the model to violate the stochastic correlation condition. But still, at the beginning, no violation of this condition occurred for a variety of values of $\rho$. Violations start appearing when ${\sigma }_{1,2}$ breaches $0.05$. Nevertheless, its rate did not surpass $9\%$ for a reasonable broad range of ${\sigma }_{1,2}$, typically up to $0.1$. Whence, this shows that condition (9) is in fact mild, at least for $ϵ$ prescribed as in (33). Also, we noticed that the violation rate decreased as we lowered the correlation parameter $\rho$.

5.2. Hand-Conducting the Model via $ϵ$

Figure 2 shows a sample path $\varrho (t,\omega )$ of the stochastic correlation coefficient as given by (5). It suggests that $\varrho$ has little variation over time. Also, from the simulations we could infer that, for any $\rho$, $99\%$ of the sample paths lay within the range of $\pm 0.10$.

We may therefore infer that our model can be a good approximation to correlated CIR models with constant correlation coefficient, with the advantage of exhibiting an analytical pricing solution.

Now, for each $\rho \in [-1,+1]$ and for $ϵ$ prescribed by the relation (33), Figure 3 exhibits the mean value over time of a corresponding sample path $\varrho$ of the stochastic correlation coefficient. Since $\varrho$ does not vary much either in time or in space, we shall rename this mean value as the operation point of the stochastic correlation coefficient, for each $\rho$.

As we can see in Figure 3, this operation point covers the interval $[-0.5,+0.5]$, which is not very broad but still a reasonable range as we vary $\rho$. Alternatively, we can depart from the maximum value of $ϵ$ given by inequality (33) to an arbitrary value in a certain confidence interval, so we can hand-conduct the curve of operation points to another more desired range.

Figure 4 shows the stochastic correlation mean surface for $ϵ$ varying in the range given by (10) and Figure 5 shows the corresponding standard deviation. These two graphs illustrate important findings. Firstly, for every $\rho$ in $(-1,1)$,

• the correlation mean vanishes when $ϵ\to 0$ (see Figure 4).
• the correlation stochasticity also vanishes when $ϵ\to 0$,(see Figure 5).

In addition, the model recast the specific case of uncorrelated CIR model when either $\rho =0$ or $ϵ=0$, for any given parameter values.

Secondly, while not near the degenerated zero value, $ϵ$ indeed introduces a correlation between the entries of the process generated by the model, whose intensity varies considerably as we hand-conduct it within the range established in (10).

These qualities underscore the model’s versatility in capturing a range of behaviors, including scenarios that are essential for theoretical analysis and practical applications.

5.3. Conditional Probability Density Functions of the Short Rate Process

We use the COS method ([5]) to recover the conditional probability density function of the $[t,T]$-integral of the short rate process R. The approximation of f is given by the following Fourier-cosine series

$$\begin{array}{c}\hfill f\left(x\right)\approx \frac{A_0}{2}+\sum _{j=1}^N{A}_{j}cos\left(j\pi \frac{x-a}{b-a}\right),x\in [a,b],\end{array}$$

(34)

for a given N. If f, with domain in $\mathbb{R}$, is a probability density function of X, then

$$\begin{array}{c}\hfill a_j\approx \frac{2}{b-a}\Re \left({exp}^{-ij\pi \frac{a}{b-a}}\widehat{f}\left(\frac{j\pi }{b-a}\right)\right)A_j,\end{array}$$

(35)

where $\widehat{f}$ is the characteristic function of X, that is

$$\begin{array}{c}\hfill \widehat{f}\left(u\right)={\int }_{\mathbb{R}}{exp}^{ixu}f\left(x\right)dx,\end{array}$$

(36)

which approximates

$$\begin{array}{c}\hfill {\int }_a^b{exp}^{ixu}f\left(x\right)dx.\end{array}$$

(37)

We set $N=100$ for the sake of truncating the series. The parameters $A_j$ are given by Equation (35) which in turn is calculated by means of the characteristic function (12). The same procedure holds for obtaining the conditional probability density function of the short rate process itself at terminal date T, in which case the characteristic function is given by Equation (21).

Setting ${\eta }_1=0$ (${\eta }_2=0$), Figure 6 (Figure 7) shows the marginal density function for $x_1\left(t\right)$ ($x_2\left(t\right)$).

Figure 8 shows the $(x_1\left(t\right),x_2\left(t\right))$-conditional probability density functions of $R\left(T\right)$ versus the correlation parameter $\rho$. High correlation, thin right tail. The lower the correlation, the more positively skewed the density is. We may also infer from the figure that the stochastic correlated CIR processes preclude negative values concerning the additive process R, so preserving the property of the single and two-factor uncorrelated square-root diffusion process (see Figure 9 for the surface of the density over time). Table 1 confirms this statement, presenting the cumulative area under the probability density function for different values of $\rho$.

Figure 10, in turn, shows the $(x_1\left(t\right),x_2\left(t\right))$-conditional probability density functions of the $[t,T]$-integral of R versus $\rho$. The lower the correlation parameter, the higher its dispersion. It is worth noticing how sensitive the shape of this density function is as we vary the parameter $\rho$, which suggests that correlation is a subject matter that cannot be neglected. The above probability density function is used for calculating the price of the IDI option.

5.4. Pricing the IDI Option

As stated by [5], the price of a European option can be approximated by

$$\begin{array}{cc}\hfill C(t,T)& \approx \frac{A_0{B}_0}{2}+\sum _{j=1}^n{A}_j{B}_j,\hfill \end{array}$$

(38)

where the $A_k$ coefficients are given by (35) and

$$\begin{array}{cc}\hfill B_j& ={\int }_a^{b}g\left(x\right)cos\left(j\pi \frac{x-a}{b-a}\right)dx,\mathrm{for}j=0,1,...,n,\hfill \end{array}$$

(39)

where $g\left(x\right)$ refers to the payoff of the financial contract.

We assume an interest rate market with underlying probability space $(\Omega ,\mathbb{F},\mathbb{Q})$ equipped with a filtration $\mathbb{F}={\left({\mathcal{F}}_t\right)}_{t\in [0,T]}$ where $\mathbb{Q}$ is the risk neutral measure.

An overnight interest rate $\overline{r}$ is the average of the interbank rate of a one-day period, calculated daily and expressed as the effective rate per annum. So, the Interbank Deposit Rate Index (IDI) accumulates discretely, according to

$$\begin{array}{c}\hfill y\left(T\right)=y\left(t\right)\prod _{j=1}^{t-1}{(1+{\overline{r}}_j)}^{\frac{1}{252}},\end{array}$$

(40)

where j denotes the end of day and ${\overline{r}}_j$ assigns the corresponding overnight rate.

If we approximate the continuously overnight rate by the instantaneous continuously compounding interest rate, i.e., $r\left(t\right)=ln(1+\overline{r}\left(t\right))$, the index can be represented by the following continuous compounding expression

$$\begin{array}{c}\hfill y\left(T\right)=y\left(t\right){exp}^{{\int }_t^{T}r\left(s\right)ds}.\end{array}$$

(41)

Given non-negative interest rates, the index is a non-decreasing function of $r\left(s\right)$.

A European call option is a contract that gives the owner the right, but not the obligation, to buy a specified amount of an underlying security at a specified price and at a specified time. The payoff of the IDI call option maturing at T is

$$\begin{array}{c}\hfill max\left(y\right(T)-K,0),\end{array}$$

(42)

where K is the strike price. Therefore, the price at time t of this option is

$$\begin{array}{ccc}\hfill C(t,T)& =& \mathbb{E}\left[{exp}^{-{\int }_t^{T}r\left(s\right)ds}max(y\left(T\right)-K,0)|r\left(t\right)\right]\hfill \\ & =& \mathbb{E}\left[max\left(y\left(t\right)-K{exp}^{-{\int }_t^{T}r\left(s\right)ds},0\right)|r\left(t\right)\right].\hfill \end{array}$$

(43)

In order to price the IDI option we compute the characteristic function of the path-dependent $[t,T]-$integral of the interest rate - namely $x(t,T)={\int }_t^{T}r\left(s\right)ds$ - which in turn will give us the values of $A_j$ associated to the model. So, the price of the IDI option can be calculated via the COS method. Equation (43) gives the payoff of the IDI call option so that the coefficients $B_j$ are given by

$$\begin{array}{c}\hfill B_0={\int }_{-ln\left(\frac{y\left(t\right)}{K}\right)}^{b}y\left(t\right)-K{\mathrm{e}}^{-x}dx=y\left(t\right)\left(ln\left({\displaystyle \frac{y\left(t\right)}{K}}\right)+b-1\right)+{\mathrm{e}}^{-b}K,\end{array}$$

(44)

and

$$\begin{array}{cc}\hfill & B_j={\int }_{-ln\left(\frac{y\left(t\right)}{K}\right)}^b\left(y\left(t\right)-K{\mathrm{e}}^{-x}\right)cos\left({\displaystyle \frac{\pi j\left(x-a\right)}{b-a}}\right)dx\hfill \\ \hfill & =\{\left(b^2-2ab+a^2\right){\mathrm{e}}^{b}y\left(t\right)sin\left(\frac{\pi jln\left(\frac{y\left(t\right)}{K}\right)+\pi aj}{b-a}\right)\hfill \\ \hfill & +(\pi a-\pi b){\mathrm{e}}^{b}jy\left(t\right)cos\left(\frac{\pi jln\left(\frac{y\left(t\right)}{K}\right)+\pi aj}{b-a}\right)\hfill \\ \hfill & +\left({\pi }^2{\mathrm{e}}^b{j}^2+\left(b^2-2ab+a^2\right){\mathrm{e}}^b\right)sin\left(\pi j\right)y\left(t\right)\hfill \\ \hfill & +\left(\left(\pi b-\pi a\right)jcos\left(\pi j\right)-{\pi }^2{j}^{2}sin\left(\pi j\right)\right)K\}\hfill \\ \hfill & .\left\{{\displaystyle \frac{(b-a){\mathrm{e}}^{-b}}{\pi j\left({\pi }^2{j}^2+b^2-2ab+a^2\right)}}\right\}.\hfill \end{array}$$

(45)

The coefficients $A_j$ are calculated with (35) by using the characteristic function given by Theorem 1.

Figure 11 gives us the IDI option prices with $T-t=1$ year and the initial value of the IDI given by $y\left(t\right)=100.000$. We note that the lower the correlation parameter, the higher the option prices, notably for at-the-money options where the low-correlated options can reach four times the price of the highly correlated contract. This is due to the effect of reducing the uncertainty implied by the probability distribution shown in Figure 10. Again, this shows that correlation is a subject that cannot be neglected.

6. Conclusion

We unveiled a bivariate model under the square-root diffusion (CIR) format equipped with a stochastic correlation coefficient for which we have obtained the analytical solution of the conditional moment generating function and the conditional characteristic function of the related diffusion process that stems from that model. Inter alia, this immediately facilitates obtaining closed-form expressions for the prices of a variety of interest rate derivatives. Remind that obtaining analytical results is a quality not achievable in the case of bivariate CIR models with constant correlation coefficient.

Varying a certain parameter ($ϵ$) serves to hand-conduct the model’s design changing the strength of correlation between the entries of the process. Figure 4 and Figure 5 illustrate this well. Also, by setting $ϵ$ equal to zero, we recover the two-factor uncorrelated square-root diffusion model.

These qualities - the analycal results in conjunction with stablishing correlations and recevering the uncorrelated CIR case - underscore the model’s versatility in capturing a range of behaviors that include scenarios that are essential for theoretical analysis and practical applications.

In addition, the simulations herein strongly suggest that the model precludes negative values, whence preserving the quality of the single and two-factor uncorrelated CIR process.

Finally, subscribed to this model, we obtained the closed-form expressions of the exact price of a zero coupon bond and a certain Asian option - the so-called IDI option.