How to handle negative interest rates in a CIR framework

In this paper, we propose a new model to address the problem of negative interest rates that preserves the analytical tractability of the original Cox-Ingersoll-Ross (CIR) model without introducing a shift to the market interest rates, because it is defined as the difference of two independent CIR processes. The strength of our model lies within the fact that it is very simple and can be calibrated to the market zero yield curve using an analytical formula. We run several numerical experiments at two different dates, once with a partially sub-zero interest rate and once with a fully negative interest rate. In both cases, we obtain good results in the sense that the model reproduces the market term structures very well. We then simulate the model using the Euler-Maruyama scheme and examine the mean, variance and distribution of the model. The latter agrees with the skewness and fat tail seen in the original CIR model. In addition, we compare the model's zero coupon prices with market prices at different future points in time. Finally, we test the market consistency of the model by evaluating swaptions with different tenors and maturities.


Introduction
The Cox-Ingersoll-Ross model (hereafter referred to as CIR model) has been regarded as the reference model in interest rate modeling by both practitioners and academics for several decades, not only because of its analytical tractability as an affine model, but also because of its derivation from a general equilibrium framework (see for example [6]), among other reasons. The well-known feature of the CIR model that ultimately led to this paper is that interest rates never become negative. This long-standing paradigm of non-negative interest rates made the CIR model and its extensions one of the most appropriate models for interest rate modeling.
Today, however, negative interest rates are very common and thus the need for models that can handle this paradigm shift is highly desirable, provided that they have as few shortcomings as possible compared to the original CIR models.
In this paper, we present a very simple and effective idea how this can be realized by modeling interest rates as the difference of two independent CIR processes, which-to the best of our knowledge-has not been considered yet.
We will propose a term structure in the risk-neutral world suitable for the difference of two independent affine processes and obtain a pricing formula for default-free zero-coupon bonds by deriving the associated Riccati equations arising from this no-arbitrage framework. In the special case of two CIR processes we will then solve the Riccati equations explicitly, which preserves the analytical tractability of its non-negative interest rate counterpart.
Afterwards, we will show some numerical experiments to demonstrate the merits of this approach in practice.
Furthermore, let the the instantaneous short-rate process be given by (1. 3) In the case where y ≡ 0, this reduces to the standard affine one-factor short rate model class.

Description of the main results
The main result consists of two main parts. First of all, we derive the zero-coupon bond price for (1.3) in the case of the difference of (1.1) and (1.2) being two independent CIR processes as in (1.4). Secondly, we provide numerical experiments to demonstrate the features of this model in Section 3. The price of a zero-coupon bond in the model r(t) = x(t) − y(t) with x and y being two independent CIR processes as in (1.4) is given by T )e By(t,T )y(t) , (1.5) where t ≤ T and for z ∈ {x, y} with φ z i ≥ 0, i = 1, 2, 3, z ∈ {x, y}, such that the Feller condition 2k z θ z ≥ σ 2 z is satisfied and The technical part of the proof is quite standard and referred to Appendix A with a description for deriving this result in Section 2. Formula (1.5) will provide the necessary ingredient for the numerical experiments in Section 3 to calibrate the model to the market term structure. In Section 3 we will perform several experiments at two different dates 30/12/2019 and 30/11/2020, where negative interest rates are observed in the market to uncover the features of this model.
There is a vast literature on interest rate modeling, among these for example the comprehensive works of [2], [5] and [11], which we cannot cover to its full extent in this small review. The most popular approach in modern interest rate modeling is the direct modeling of short rates r(t) under a risk-neutral measure Q inspired by no-arbitrage arguments, especially, because the price at time t > 0 of a contingent claim with payoff H T , T > t, under the risk-neutral measure is given by (cf. [22]) where E Q t denotes the conditional expectation with respect to some filtration F t under measure Q. In particular, choosing H T := P (T, T ) = 1, where P (t, T ) denotes a zero-coupon bond, gives rise to a convenient way to calibrate a short rate model to the market term structure, which we will utilize for our approach as well.
Starting with the pioneering works of Merton [17] in 1973 and Vasicek [23] in 1977, many one-factor short rate models were introduced, see [5] for a detailed overview. Among all, a model that has had a particular importance in the past, equally among both practitioners and academics, is the well-known CIR model, proposed by Cox, Ingersoll & Ross in [7]. It provides the basis for this paper and is a generalization to the Vasicek model by introducing a non-constant volatility given by (1.4).
Clearly, the square root term precludes the possibility of negative interest rates and under the assumption of the Feller-condition 2kθ ≥ σ 2 , see for instance [14], the origin is inaccessible.
These two properties combined with its analytical tractability make the CIR model well-suited for a non-negative interest rate setting.
There is a rich literature on extensions to the classical CIR model in order to obtain more sophisticated models, which could fit the market data better, allowing to price interest rate derivatives more accurately. For example, Chen in [1] proposed a three-factor model; Brigo and Mercurio in [5] proposed a jump diffusion model (JCIR). In order to include time dependent coefficients in (1.4), Brigo and Mercurio in [4] proposed to add a deterministic function into equation (1.4). This model, called CIR++, is able to fit the observed term structure of interest rates exactly, while preserving the positivity of the process r(t). Brigo and El-Bachir in [3] generalized the CIR++ model by adding a jump term described by a time-homogeneous Poisson process and Brigo and Mercurio in [5] studied the CIR2++ model. Another way to generalize the CIR model by including time dependent coefficients in equation (1.4) was introduced by Jamshidian in [13] and Maghsoodi in [15], which are known as extended CIR models.
But in the last decade the financial industry encountered a paradigm shift by allowing the possibility of negative interest rates, making the classical CIR model unsuitable.
One way to handle the challenges entailed by negative interest rates is to use Gaussian models with one or more factors, such as the Hull and White model (see [12]), which also has a very good analytical tractability. A generalization of these models with a good calibration to swaption market prices was found in [9], while Mercurio and Pallavicini in [16] proposed a mixing Gaussian model coupled with parameter uncertainty.
But the glamor of the CIR model is still alive even in the current market environment with negative interest rates. Orlando et al. suggest in several papers (cf. [19], [20] and [21]) a new framework, which they call CIR# model, that fits the term structure of interest rates. Additionally, it preserves the market volatility, as well as the analytical tractability of the original CIR model. Their new methodology consists in partitioning the entire available market data sample, which usually consists of a mixture of probability distributions of the same type. They use a technique to detect suitable sub-samples with normal or gamma distributions. In a next step, they calibrate the CIR parameters to shifted market interest rates, such that the interest rates are positive, and use a Monte Carlo scheme to simulate the expected value of interest rates.
In this paper, however, we introduce a new methodology for handling the challenges arising from negative interest rates. In our model, the instantaneous spot rate is defined as the difference between two independent classical CIR processes, which allows the preservation of the analytical tractability of the original CIR model without introducing any shift to the market interest rates.
The paper is organized as follows. In Section 2 we introduce the model in a general affine model setup and describe our main result Theorem 1.1. We will derive the Riccati equations associated with the proposed term structure suitable for the difference of two independent affine processes and solve those explicitly in a CIR framework.
After that, in Section 3, we will conduct some numerical experiments. First, we calibrate our model via (1.5) to the market data at 30/12/2019 and 30/11/2020 in Section 3.2. Subsequently, we simulate the model by using the Euler-Maruyama scheme in Section 3.3 and study the mean, variance and distribution of the model in Section 3.4. Then we test how the calibrated model performs when pricing zero-coupon bonds at future times in Section 3.5 and conclude our numerical tests by pricing swaptions in Section 3.6. Finally, we summarize the results of the paper in Section 4 and discuss possible extensions for future research.

A model for negative interest rates
We will now describe how Theorem 1.1 can be derived. As aforementioned, we consider all dynamics under the risk-neutral measure Q and give now a heuristic argument, why it makes sense to choose the term structure in Theorem 1.1 as in (1.5).
Suppose, that x(t) and y(t) are both independent affine processes. Then the the price of a zero-coupon bond for each of them separately (cf. [5] p. 69) is given by where z ∈ {x, y} and E Q t denotes the conditional expectation with respect to F t under the measure Q. Now, consider r(t) = x(t) − y(t), then we have by linearity and independence If we concentrate in (2.1) only on the right-hand side, it would make sense for two independent processes x and y that we can apply these formulas with a change of sign in front of B y , leading to In the following Lemma we will make this argument rigorous.

Then, the price of a Zero-coupon bond is given by
where A z and B z , z ∈ {x, y}, are deterministic functions and are a classical solution to the following system of Riccati equations Introducing dependence between x and y suggests a coupling of A x and B x to A y and B y and might have an impact on the analytical tractability, but is left for future research.
It is well-known that, under the conditions 2k z θ z ≥ σ 2 z , z ∈ {x, y}, the processes x(t) and y(t) are strictly positive for every t ≥ 0 (see for instance [7] or [14]). We underline that even if the processes x(t) and y(t) are strictly positive, the instantaneous spot rate r(t) could be negative since it is defined as the difference of x(t) and y(t) for every t > 0, which is illustrated in Figure 1 together with several percentiles of r(t).  Table 2.

Numerical tests
We will now perform some numerical experiments in our model. In Section 3.1 we will briefly discuss the market data, which we will use to perform all numerical tests in the subsequent sections. Afterwards, we will describe the calibration procedure of our model to the zerocoupon curves at 30/12/2019 and 30/11/2020 in Section 3.2. This is followed by a short subsection on simulating the model with the Euler-Maruyama scheme in Section 3.3 and we investigate the mean, variance and distribution of the short rate model in Section 3.4. In Section 3.5 we price zero-coupon bonds at future dates and compare the results to the market prices. Last but not least, Section 3.6 will show results on pricing swaptions in our model.
We used for the calculations Matlab 2021a with the (Global) Optimization Toolbox running on Windows 10 Pro, on a machine with the following specifications: processor Intel(R) Core(TM) i7-8750H CPU @ 2.20 GHz and 2x32 GB (Dual Channel) Samsung SODIMM DDR4 RAM @ 2667 MHz.

Market Data
To obtain the market zero-coupon bond term structure, we first build the EUR Euribor-swap curve which is created from the most liquid interest rate instruments available in the market and constructed as follows: We consider deposit rates and Euribor rates with maturity from one day to one year and par-swap rates versus six-month Euribor rates with maturity from two years to thirty years. Then the zero interest curve and the zero-coupon bond curve are calculated using a standard "bootstrapping" technique in conjunction with cubic spline interpolation of the continuously compounded rate (cf. [18] for more details).
We choose two different dates and we take the data at the end of each business day. In particular, we test our model at 30/12/2019 and at 30/11/2020. At the first date, the zero interest rates were negative up to year six, while at the second date the entire zero interest rate structure was negative. In Table 6 and in Table 7 we report the zero interest rate curve and the zero-coupon bond curve at the two different dates.
Furthermore, in Section 3.6 we need the strikes to compute the model swaption prices and the market prices of the swaptions to compare our model, which are for both dates in the Appendix in Table 8, Table 9, Table 10 and Table 11, respectively. All data has been downloaded from Bloomberg and is used in the following subsections for our numerical experiments. We start in the next subsection with calibrating our model to the zero-coupon curve.

Calibration
In this subsection we will discuss how we calibrate our model to the market zero-coupon curve given in Table 6 and Table 7 by using the formula derived in (1.5).
Let us denote Π : We will formulate the calibration procedure as a constraint minimization problem in R 8 for the parameters Π with objective where n ∈ N is the number of time points, where market data is available, and T i , i = 1, . . . , n are these maturities. The market zero-coupon curve is denoted by P M (0, T i ) and P (Π; 0, T i ) is the price of a zero-coupon bond in our model given by (1.5) with parameters Π.
The objective function describes the relative square difference between the market zerocoupon bond prices and the theoretical prices from the model given by (1.5).
The set of admissible parameters A will consist of the following constraints arising from the well-definedness of the formulas (1.7): (i) First of all, let us note that there is a one-to-one correspondence between the parameters Π and k z , σ z and θ z if one is looking for positive real solutions only. We have (v) A positive mean for each CIR process, i.e. θ z ≥ 0, is by positivity of σ 2 z and k z equivalent to φ z 3 ≥ 0, which is already satisfied by the Feller condition; (vi) The parameter φ z 1 , assuming that it is real-valued, is positive by definition, meaning that by the positivity of the mean reversion speed, φ z 2 will be as well. Therefore, all φ are positive; (vii) As both CIR processes x t and y t , individually, are positive processes, we additionally require x 0 ≥ 0 and y 0 ≥ 0.
The advantage of using the parameters Π instead of k z , σ z and θ z is that we can rewrite these conditions as a system of linear inequality constraints in matrix notation A · Π ≤ 0, where In total, the set of admissible parameters is given by Finally, a solution Π * to the calibration problem is a minimizer of To solve (3.4) numerically, we want to use Matlab's function fmincon in the (Global) Optimization Toolbox. In order to use this function, we need an initial guess of the parameter Π and the computational time will depend on that choice. In Table 1 we present a few choices for initial guesses of Π. The first row for each date 30/12/2019 or 30/11/2020 refers to Matlab's function ga in the (Global) Optimization Toolbox, which uses a generic global optimization algorithm to find a solution of (3.4) without starting from an initial guess, which takes a long time to compute, roughly 35 to 43 seconds. In the following three rows are three manual initial guesses. We can see that the first two choices work for both dates exceptionally fast (0.3 seconds) and the accuracy is almost identical to all other choices, making this model a good choice if live calibration to the data is needed, which we also use in the following numerical experiments. In the last two rows we used random starting parameters to demonstrate that the error remains stable but the computational time varies.
For the algorithms used by Matlab we refer the reader to [10], in the context of financial mathematics.
The results of the aforementioned calibration procedure are displayed in Table 2

Euler-Monte-Carlo simulation
In order to forecast the future expected interest rate, we use the Euler-Maruyama scheme to simulate the instantaneous spot rate r (1.3). We refer to [8] and the references therein for a list of different Euler-type methods to simulate a CIR process. In our experiments, we simulate the processes x(t) and y(t) by the truncated Euler scheme defined as follows: First of all, we fix a homogeneous time grid 0 = t 0 ≤ t 1 ≤ · · · ≤ t N = T for the interval [0, T ] with N + 1 time points and mesh ∆t i := t i+1 − t i ≡ ∆ := T N for all i = 0, . . . , N − 1. Secondly, we simulate the two independent Brownian motions W z , z ∈ {x, y}, and define their time We choose the max inside the square-root to ensure that the square-root remains real, because due to discretization effects the positivity of x(t i ) and y(t i ) might be violated.    (Table 6) to the mean of the model discount factors with absolute errors at 30/12/2019 with parameters given in Table 2  difference of the mean over all simulations of our model to the market data, at each maturity can be found in the appendix in Table 5.

Mean and variance
The F s -conditional mean and variance of the CIR process are well-known (cf. [5] p. 66 equation (3.23)) and are given by where z ∈ {x, y}. In the case of the difference of two CIR processes we have and by independence In Figure 5 we show for each date the histogram of the short rate distribution after 30 years.
To describe the distribution of r(t) after 30 years better, we also compare it to the density of a normal random variable with the same mean and variance. As one expects, the distribution of r shows a slight skewness and fatter tail with respect to the normal distribution.   Table 2, obtained with (3.6) and (3.7). The left picture shows the results at 30/12/2019 and the right at 30/11/2020.

t-Forward zero-coupon prices
In this subsection we examine, if our model is able to replicate the interest rate term structure not only at the start date but also in the future. Therefore, we compare the market and model prices of zero-coupon bonds at future times t by using the continuously-compounded spot interest rates (cf. [5] Definition 1.2.3.) to obtain the market prices.
To be more precise, let us briefly recall the definition of the continuously-compounded spot interest rates where we use the day-count convention T − t in years.
To compute this rate for the market at a fixed future date t > 0 with maturity T > t we use where R M (0, t) and R M (0, T ) are the market zero rates with maturity t and T respectively.
By rearranging (3.8) we derive the market zero-coupon price at a future date t with maturity   Figure 6.: Mean and 99 % confidence interval of the model t-forward zero-coupon prices compared to the market t-forward zero-coupon prices at t = 1 using the calibrated parameters in Table 2  In Figure 6, Figure 7 and Figure 8 we show a comparison of (3.10) to (1.5) at 30/12/2019 and 30/11/2020 with parameters Π * from Table 2 for future times t = 1, 3, 5, respectively.
The behavior shown by the model t-forward prices coincides with typical one-factor short interest rate models. At short future dates, one year for example, the model is able to reproduce the market forward zero-coupon price for 30/12/2019, i.e. the error is of magnitude 0.01, whereas at long future dates, e.g. Figure 8, the model prices are deviating further from the market prices, i.e. the error is of magnitude 0.03. For 30/11/2020 we can see that the errors are not very large but the model seems to have difficulties to match the shape of the predicted market prices.

Pricing swaptions
In this subsection we test if our model is market consistent, in the sense whether the model is able to reproduce market swaption prices or not.
We compare market swaption prices to model swaption prices with different tenors (1, 2, 5, 7, 10 years) and maturities (1, 2, 5, 7, 10, 15, 20 years Figure 8.: Mean and 99 % confidence interval of the model t-forward zero-coupon prices compared to the market t-forward zero-coupon prices at t = 5 using the calibrated parameters in Table 2, ∆ = 1 256 and M = 10000. The left picture shows the results at 30/12/2019 and the right at 30/11/2020.  Table 3 and Table 4, respectively. We notice that, similar to one-factor short interest rate model, our model fails to capture the full swaption volatility surface. This result is not surprising, since the model uses essentially a single volatility factor due to the fact that the model parameters are constant and the Brownian motion are independent. A way to make the model able to match the entire volatility surface is the inclusion of time-dependent parameters to fit exactly the market term structure or, equivalently, adding a deterministic function of time to the short rate process r, which is left for future research. That being said, the authors would like to stress that this is a first step in this methodology of using two CIR processes as a difference to model negative interest rates. Naturally, its extensions, such as considering time-dependent coefficients to fit the market term structure perfectly or, equivalently, adding a deterministic shift extension in the sense of [5], will be a next step and is left for future research.

A. Derivation of the Riccati equations and coefficients
Let everything be as in Lemma 2.1. In particular, let To derive the Riccati equations (2.3) we use the fact, that we are modelling under the martingale measure Q, therefore the discounted price process exp − t 0 r(s)ds P (t, T ) needs to be a martingale.
By independence of x and y, as well as Itô's formula we derive after some algebra Now, in order to be a martingale the parts of bounded variation have to vanish, which leads us after rearranging the terms to Thus, we derive the following Riccati System A solution to the last equation can be found by further assuming that the individual x and y parts will be zero, leading to two separate equations the same way, leading to 2 k 2 − 2 h k + 2 k 2 e 2 h τ − 2 σ 2 − 4 σ 2 e h τ − 2 σ 2 e 2 h τ + 2 h k e 2 h τ .
In total, we see that the denominator differs only by a sign, hence ∂ t B + 1 2 σ 2 B 2 − kB = −1, which yields the claim.
Verification for A The formula can be derived by just integrating and taking the exponential.
Let us just take the logarithm and derivative to verify this formula: Now, taking the derivative yields After undoing the substitution for τ this is equal to kθB, which yields the claim.
By (1.5) we therefore have

Conflicts of interests
The authors have no relevant financial or non-financial interests to disclose.

Data availability
All data generated or analysed during this study are included in this published article. In particular the code to produce the numerical experiments is available at https://github.com/kevinkamm/CIR-.