Welcome to THETAWIKI. If you like to create or edit a page please make sure to login or register an account. All registered users please make sure to provide a valid email address.

# CIR model

## The CIR model

The Cox-Ingersoll-Ross (CIR) (1985) model[1] describes the evolution of short rates as having the following dynamics:

$dr(t) = a\, \left [ \theta \,-\, r(t) \right ] dt \,+\, \sigma \sqrt{r(t)} \, dW(t), \qquad (CIR.1)$

Where the terms $a, \, \theta , \, \sigma$ are positive constants, and denote respectively, the mean reversion speed, the mean reversion level and volatility of short rate $r(t)$. For our convenience, we assume the above short rate process is under the risk-neutral measure $Q$, i.e. the term $W(t)$ is a standard $Q$ wiener process. We can, however, easily change the above process from the risk-neutral measure $Q$ to another measure, such as the physical measure $P$, defined on the same sample space.

The CIR (1985) model is an equilibrium model, in that the initial term structure is an output rather than an input. As such, we need to calibrate to the current zero coupon bond prices first before pricing bond options. However, exact fitting the above model to an array of zero coupon bond prices may not be possible, as there are two few numbers of parameters for more numbers of bond prices. The advantage of having the CIR (1985) model is that analytical formulas for zero coupon bond price and zero coupon bond option prices are available, and that short rates are guaranteed to be positive.

To obtain future short rate trajectories, we discretize in time equation (CIR.1) before using various numerical schemes for simulations. For our sample project, we choose to use the following simulation schemes:

• modified implicit Milstein scheme
• Andersen(2008)’s quadratic exponential (QE) scheme

In the following, we will first document the above simulation schemes, then show how we implement them in ThetaML.

## The Modified implicit Milstein scheme for the CIR model

We time discretize equation (CIR.1) on the time interval $\left [0, \, T \right ]$ with $n$ equal time steps on the set of time schedules $0 \, = \, t_0 \, \leq \, t_1 \, \leq \, \cdot \, \leq \, t_n \, = \, T$.

The explicit Euler scheme of equation (CIR.1) with the above time discretization is given as follows:

$\hat{r}(t \,+ \, \Delta) = \hat{r}(t) \, + \, a\, \left [ \theta \,-\,\hat{ r}(t) \right ] \Delta\,+\, \sigma \sqrt{\hat{r}(t)} \, \sqrt{\Delta} \,Z(t), \qquad (CIR.2)$

where $\hat{r}(t)$ is the discretized short rate evaluated at the left hand side of discretization time interval, and the term $Z(t)$ is a Gaussian random variable. The above explicit Euler scheme, however, is not complete, since negative short rate can not go into the $\sqrt{ \, }$ sign. We therefore modify the discretized equation (CIR.2) as follows:

$\hat{r}(t \,+ \, \Delta) = \hat{r}(t) \, + \, a\, \left [ \theta \,-\,\hat{ r}(t)^+ \right ] \Delta\,+\, \sigma \sqrt{\hat{r}(t)^+} \, \sqrt{\Delta} \,Z(t), \qquad (CIR.3)$

as such, negative short rate $\hat{r}(t)$ is possible. Typical of interest rate market, the Feller condition $2a \theta \, > \, \sigma^2$ is satisfied. The explicit Euler scheme in equation (CIR.3) has first-order weak convergence, to improve the convergence rate, we use the implicit Milstein scheme, as what documented below.

The implicit Milstein scheme using the above discretization time schedule is as follows:

$\hat{r}(t \,+ \, \Delta) = \hat{r}(t) \, + \, a\, \left [ \theta \,-\,\hat{r}(t \,+ \, \Delta) \right ] \Delta\,+\, \sigma \sqrt{\hat{r}(t)} \, \sqrt{\Delta} \,Z(t),$
$\hat{r}(t \,+ \, \Delta) \left ( 1\,+ \,a \, \Delta \right ) = \hat{r}(t) \, + \, a\, \theta \,\Delta\,+\, \sigma \sqrt{\hat{r}(t)} \, \sqrt{\Delta} \,Z(t),$
$\hat{r}(t \,+ \, \Delta) =\frac{\hat{r}(t) \, + \, a\, \theta \,\Delta\,+\, \sigma \sqrt{\hat{r}(t)} \, \sqrt{\Delta} \,Z(t)}{ \left ( 1\,+ \,a \, \Delta \right ) }. \qquad (CIR.4)$

The implicit Milstein scheme in equation (CIR.4) has strictly positive rates if the condition $4a \theta \, > \, \sigma^2$ holds, in cases where this condition is not satisfied, we turn to equation (CIR.3) for updating short rate process.

## Andersen (2008)’s quadratic exponential (QE) scheme for the CIR model

Andersen (2008)’s QE scheme uses a combination of squared Gaussian and exponential distributions.

Case 1.

When $r(t)$ is large, the QE scheme moment matches the first two moments of CIR short rate process to the first two moments of the chi-square distribution.

$\hat{r}(t \,+ \, \Delta) = a \, \left ( b \,+\, Z_r \right ) ^2, \qquad (CIR.5)$

where the parameters $a, \, b$ are to be determined from the moment matching process, and $Z$ is a standard Gaussian random variable. The above formulation comes from the observation that a squared Gaussian random variable is $\chi^2$ distributed.

We now turn to the moment matching scheme to find a representation for the parameters $a, \, b$, respectively.

To obtain the first two conditional moments of the CIR (1985) model, we solve for the integral representation of equation (CIR.1) from $t$ to $t\,+\, \Delta$ as follows.

$dr(t)\,+\,a\,r(t)\, dt = a\, \theta \, dt \,+\, \sigma \sqrt{r(t)} \, dW(t)$
$d\left (e^{at} \, r(t) \right ) = a\, \theta \, e^{at} \,dt \,+\, \sigma \, e^{at} \, \sqrt{r(t)} \, dW(t)$
$e^{a\, (t\,+\, \Delta)} \, r(t\,+\, \Delta) \,-\, e^{at} \, r(t) = a\, \theta \, \int_t^{ t\,+\, \Delta} e^{au} \,du \,+\, \sigma \, \int_t^ {t\,+\, \Delta} e^{au}\, \sqrt{r(u)} \, dW(u)$
$r(t\,+\, \Delta) = r(t) e^{-a\, \Delta} \, + \, \theta \, \left ( 1\,-\,e^{-a\, \Delta} \right )\,+\, \sigma \, \int_t^{ t\,+\, \Delta} e^{-a\,( t\,+\, \Delta -\,u)}\, \sqrt{r(u)} \, dW(u) . \qquad (CIR.6)$

From equation (CIR.6), we have the first two conditional moments for the CIR(1985) short rate model:

$\mathbb{E} \left [ r(t\,+\, \Delta ) \mid r(t) \right ] = r(t)e^{-a\, \Delta } \,+\, \theta \, \left ( 1\,-\,e^{-a\, \Delta} \right ) := m, \qquad (CIR.7)$
$Var \left [ r(t\,+\, \Delta ) \mid r(t) \right ] = \frac{\sigma^2}{a} r(t)e^{-a\, \Delta } \left ( 1\,-\,e^{-a\, \Delta} \right ) \,+\, \frac{\sigma^2}{2a} \theta \, \left ( 1\,-\,e^{-a\, \Delta } \right )^2 := s^2. \qquad (CIR.8)$

The first two moments for equation (CIR.5) is:

$\mathbb{E} \left [r(t\,+\, \Delta) \mid r(t) \right ] = a\, \left (1\,+\,b^2 \right ), \qquad (CIR.9)$
$Var \left [ r( t\,+\, \Delta ) \mid r(t) \right ] = 2\, a^2 \, \left ( 1\,+\,2\,b^2 \right ). \qquad (CIR.10)$

Equating the right hand sides of equations (CIR.7) and (CIR.8) to those of equations (CIR.9) and (CIR.10) respectively, we have

$a\, \left ( 1\,+\,b^2 \right ) = m, \qquad (CIR.11)$
$2\, a^2 \, \left ( 1\,+\,2\,b^2 \right ) = s^2. \qquad (CIR.12)$

Let $\phi ^{-1} = \frac{m^2}{s^2},$ and plug in equations (CIR.11) and (CIR.12), we have

$\phi ^{-1} = \frac{m^2}{s^2},$
$= \frac{ a^2\, \left ( 1\,+\,b^2 \right ) ^2 }{2\, a^2 \, \left ( 1\,+\,2\,b^2 \right )},$

Denoting $x = b^2$, we arrive at the following equation

$\phi ^{-1} = \frac{ a^2\, \left ( 1\,+\,b^2 \right ) ^2 }{2\, a^2 \, \left ( 1\,+\,2\,b^2 \right )},$
$\Rightarrow x^2 \,+\, 2 \left ( 1 \,-\, 2 \phi ^{-1} \right ) \, x \,+\, 1\,-\, 2\, \phi ^{-1} = 0. \qquad (CIR.13)$

Solving equation (CIR.13) for $x$, we have

$b^2 = x = - \left ( 1 \,-\, 2 \phi ^{-1} \right ) \pm \sqrt{ 2\, \phi ^{-1} \, \left (2\, \phi ^{-1} \,-\, 1 \right ) }. \qquad (CIR.14)$

Having solved for $b^2$, the expression for $a$ can be easily obtained from equation (CIR.11) as

$a = \frac{m}{1\,+\,b^2}. \qquad (CIR.15)$

Case 2.

When $r(t)$ is low, Andersen (2008) proposes using the following approximation for the density of $\hat{r}(t\,+\, \Delta)$:

$Pr \left ( \hat{r}(t\,+\, \Delta) \, \in \, \left [r, \, r\,+\,dr \right ] \right ) \, \approx \, \left ( p\, \delta{(r)} \,+\, \beta \, (1\,-\,p)\,e^{-\beta \, r} \right ) \, dr, \, r \, \geqq \, 0, \qquad (CIR.16)$

Where $\delta{(r)}$ is the Dirac delta function, and the terms $p, \, \beta$ are non-negative constants to be determined.

The cumulative distribution of equation (CIR.16) is

$\Phi \left ( r \right ) = Pr \left ( \hat{r}(t\,+\, \Delta) \, \leqq r \right ) = p\, + \, (1\,-\,p)\, \left ( 1 \,-\,e^{-\beta \, r} \right ). \qquad (CIR.17)$

Let $\phi, \, m,\, s$ be as defined as in Case 1, Andersen (2008) suggests setting $p$ as follows

$p = \frac{ \phi \,-\, 1}{\phi \,+\, 1} \in \left [0, \, 1 \right ], \qquad (CIR.18)$

and

$\beta = \frac{ 1 \,-\, p}{m } = \frac{2}{m\, \left ( \phi \,+\,1 \right ) } \, > \, 0. \qquad (CIR.19)$

The QE algorithm can be found in Andersen (2008)[2].

In the following, we give code examples for the modified implicit Milstein scheme and Andersen (2008)’s QE scheme[2].

## ThetaML implementation of the modified implicit Milstein scheme of the CIR model

This ThetaML code example simulates CIR (1985) short rates using the modified implicit Milstein scheme.

model CIR_Milstein
%this model uses the implicit Milstein scheme to simulate the square-root CIR process
%Feller condition: 2*a*tht > sigma^2, then the process for short rate can never reach zero,
%else, the origin is accessible.
import a      "Mean reversion speed"
import tht    "Mean reversion level"
import sigma  "CIR short rate volatility"
import r0     "Initial short rate"
import ht     "Discretization time step"
import t      "Simulation time horizon"
export rt     "Simulated short rates using the implicit Milstein scheme"

rt = r0           %initialize short rate
loop t/ht
theta ht        %the 'theta' command passes time by 'ht' time units
Z = randn()
if (4*a*tht > sigma^2)
%if this condition holds, we have strictly positive short rate
rt = 1/(1+a*ht) * ( rt + a*tht*ht + sigma * sqrt(rt) * sqrt(ht) * Z
+ 1/4 * sigma^2 * ht*(Z^2 - 1) )
else
%else, short rate can be negative, use the following modified equation
rt = rt + a * (tht - max(rt, 0))*ht + sigma * sqrt(max(rt,0)) * sqrt(ht)*Z;
end
end
end

The following graph shows simulated CIR (1985) short rates using the modified implicit Milstein scheme. The model parameters are as follows:

• $a = 0.5$, “Mean reversion speed”
• $\theta = 0.03$, “Mean reversion level”
• $\sigma = 0.05$, “CIR short rate volatility”
• $ht = 1/12$, “Discretization time step”
• $t = 10$, “Simulation time horizon”
• $n = 10000$, “Number of simulations”

This graph shows the histogram of simulated CIR (1985) short rates at time $T = 10$, using the modified implicit Milstein scheme.

## ThetaML implementation of Andersen (2008)’s QE scheme of the CIR model

This ThetaML code example simulates CIR (1985) short rates using Andersen (2008)’s QE scheme[2].

model CIR_QEscheme
%this model uses Andersen's QE scheme to simulate square root short rate process
%Andersen's QE scheme uses moment matching scheme to match the first two conditional moments
%of the CIR type process to the first two moments of chi-square distribution
import a      "Mean reversion speed"
import tht    "Mean reversion level"
import sigma  "CIR short rate volatility"
import r0     "Initial short rate"
import ht     "Discretization time step"
import t      "Simulation time horizon"
export rt     "Simulated short rates using Andersens QE scheme"

rt = r0       %initialize the short rate
%constant for the switching rule used in the QE scheme
shi_crt = 1.5
loop t/ht
theta ht
%conditional mean of the CIR model
m = tht + (rt-tht) * exp(-a*ht)
%conditional variance of the CIR model
s2 = rt * sigma^2/a * exp(-a*ht) * (1-exp(-a*ht)) + tht*sigma^2/(2*a) * (1-exp(-a*ht))^2
%switching rule, Andersen's paper, equation (19)
shi = s2/(m^2)

%Andersen's QE algorithm 3.2.4, p16 - 17
u = rand()
Z1 = norminv(u,0,1)
if shi <= shi_crt
%for sufficiently large value s of rt
c4 = 2/shi
b2 = max( c4 - 1 + sqrt(c4*(c4-1)), 0 )
aa = m/(1+b2)
rt = aa * (sqrt(b2)+Z1)^2
else
%for low values of rt
p = (shi-1)/(shi+1)
beta = (1-p)/m
if (u <= p)
rt = 0
else
rt = log((1-p)/(1-u)) / beta
end
end
end
end

This graph shows simulated CIR (1985) short rates using Andersen (2008)’s QE scheme[2]. The model parameters are as follows:

• $a = 0.5$, “Mean reversion speed”
• $\theta = 0.03$, “Mean reversion level”
• $\sigma = 0.05$, “CIR short rate volatility”
• $ht = 1/12$, “Discretization time step”
• $t = 10$, “Simulation time horizon”
• $n = 10000$, “Number of simulations”

This graph shows the histogram of simulated CIR (1985) short rates at time $T = 10$, using Andersen (2008)’s QE scheme.

## References

1. Cox, J.C., Ingersoll, J.E. and Ross, S.A., 1985, A Theory of the Term Structure of Interest Rates, Econometrica 53, p385-407.
2. Andersen, L., 2008, Simple and Efficient Simulation of the Heston Stochastic Volatility Model, Journal of Computational Finance, 11 (3), p1-42.