Approximation (Milstein Scheme) of the Cox-Ingersoll-Ross process

January 16, 2023

The CIR process is given by the SDE:

(1) $\begin{equation*} d r(t)=\alpha(b-r(t)) d t+\sigma \sqrt{r(t)} d W(t) \end{equation*}$

(cf. Glasserman, p. 108); where W(t) is a Brownian motion. Take α=0.2, σ=0.02, b=0.05, and r(0) = 0.044.

By application of the first-order Milstein scheme to the CIR model, we have:

(2) $\begin{equation*} \begin{aligned} r\left(t_{i+1}\right)= & r\left(t_i\right)+\alpha \left(b-r\left(t_i\right)\right)\left(t_{i+1}-t_i\right) \\ & +\sigma \sqrt{r(t_i)} \sqrt{t_{i+1}-t_i} Z_{i+1}+\frac{1}{4} \sigma^2\left(t_{i+1}-t_i\right)\left(Z_{i+1}^2-1\right) \end{aligned} \end{equation*}$

where Z~N(0, 1), so the Monte Carlo simulation could be:

f <- function(non_par){
    alpha = 0.2
    sigma = 0.02
    b = 0.05
    r0 = 0.044
    q0 = 1
    s0 = 0
    m = 1000
    dt = 5 / m
    p = 1
    
    Z = rnorm(m, mean = 0, sd = 1)
    
    r <- as.numeric(list())
    r[1] = r0 + alpha*(b-r0)*dt + sigma*sqrt(r0)*sqrt(dt)*Z[1] + 1/4*sigma**2*dt*(Z[1]**2 - 1)
    for(i in 2:m){
        r[i] = r[i-1] + alpha*(b-r[i-1])*dt + sigma*sqrt(r[i-1])*sqrt(dt)*Z[i] + 1/4*sigma**2*dt*(Z[i]**2 - 1)
    }
    
    q <- as.numeric(list())
    q[1] = q0 + q0*r0*dt
    for(i in 2:m){
        q[i] = q[i-1] + q[i-1]*r[i-1]*dt
    }
    
    s <- as.numeric(list())
    s[1] = s0 + (p + s0*r0)*dt
    for(i in 2:m){
        s[i] = s[i-1] + (p + s[i-1]*r[i-1])*dt
    }
    return(c(s[m], q[m]))
}

num <- 1000
test <- c(1:num)
Y <- data.frame(lapply(test, f));
sum_s = 0
sum_q = 0
for (i in 1:num){
    sum_s <- sum_s + Y[[i]][1]
    sum_q <- sum_q + Y[[i]][2]
}
p <- sum_q / sum_s
p

The simulation part is stored in {r(t)}, while other parts are used for (Asmussen and Glynn, p. 291, Problem 4.1).

Tags:Monte Carlo