On weak approximations of CIR equation with high volatility

Abstract

We propose two new positive weak second-order approximations for the CIR equation $\mathrm{d}{X}_t=(a-b{X}_t)\mathrm{d}t+\sigma \sqrt{X_t}\mathrm{d}{B}_t$ based on splitting, at each step, the equation into the deterministic part $\mathrm{d}{X}_t=(a-b{X}_t)\mathrm{d}t$, which is solved exactly, and the stochastic part $\mathrm{d}{X}_t=\sigma \sqrt{X_t}\mathrm{d}{B}_t$, which is approximated in distribution. The schemes are illustrated by encouraging simulation results.

Introduction

In this paper, we are interested in positive weak second-order approximations of the solution of CIR equation (Cox–Ingersoll–Ross [5])

$$\mathrm{d}{X}_t=(a-b{X}_t)\mathrm{d}t+\sigma \sqrt{X_t}\mathrm{d}{B}_t\text{,}{X}_0=x_0\ge 0\text{,}$$

where $B$ is a standard Brownian motion, and a, $\sigma \ge 0$, $b\in R$. Square-root diffusions arise in applications ranging from finance to medicine; see, e.g., Gillespie [7] (chemical kinetics), Fox [6] (neuroscience), Allen [4, Chap. 5] (population dynamics), among others. However, the CIR equation, which was initially introduced to model the short interest rate [5], seems to be most widely used in finance and financial mathematics. Mathematically, its main qualitative features are the positivity of the solution and known analytic expressions for its moments. Note that the distribution (of increments) of CIR process is known explicitly as a noncentral ${\chi }^2$ distribution (see, e.g., Lamberton and Lapeyre [9, pp. 130–131]), and thus its exact simulation is possible; however, it is too slow in comparison with discretization schemes. Therefore, construction and analysis of positivity-preserving approximation methods for CIR and other square-root diffusion equations is still a topic of great interest; see, e.g., Halley et al. [8].

The main problem with developing numerical methods for square-root diffusions is $\dots$ the square-root itself, which has unbounded derivatives near zero. Therefore, discretization schemes that (explicitly or implicitly) involve the derivatives of the coefficients – even when they assure the positivity of approximation – usually loose their accuracy near zero, especially, for large $\sigma$; the larger $\sigma$ is, the more concentrated near zero the value-distributions of CIR process are. One way to get round this difficulty is modifying the scheme considered by switching near zero to another scheme, which (1) is sufficiently “smooth” and (2) sufficiently accurate near zero; we refer to Alfonsi [3] and references therein.

Recall that a family of processes $\{X^h\text{,}h>0\}$ is said to be an approximation of the solution $X$ in the weak sense (or, shortly, a weak approximation of $X$) of order $n$ on the time interval $[0\text{,}T]$ if

$$\text{E}f(X_T^h)-\text{E}f(X_T)=\text{O}(h^n)\text{,}h\to 0\text{,}$$

for a rather wide class of (test) functions $f$. We consider the approximations of the form

$$X_0^h=x\text{,}{X}_{(k+1)h}^h=A(X_{kh}^h\text{,}h\text{,}\Delta B_k)\text{,}k=0\text{,}1\text{,}2\text{,}\dots \text{,}$$

where the (increment) function $A(x\text{,}h\text{,}y)\text{,}(x\text{,}h\text{,}y)\in (0\text{,}\infty )\times [0\text{,}\infty )\times IR$, is such that $A(x\text{,}0\text{,}0)=x$ and $\Delta B_k=B_{(k+1)h}-B_{kh}$. We suppose them to be extended from the grid $\{kh\text{,}k=0\text{,}1\text{,}\dots \text{,}\}$ the whole interval $[0\text{,}T]$ as a step function, i.e., $X_t^h:=X_{[t{h}^{-1}]h}^h$, where $[\cdot ]$ denotes the integer part. To show the starting point $X_0^h=x$, we shall write $X^h(x\text{,}t)$ instead of $X_t^h$. Thus, in our setting, $X^h(x\text{,}h)=A(x\text{,}h\text{,}{B}_h)$; therefore, with some abuse of language, we identify increment functions with the approximations they define. Alfonsi [1] constructed and studied several first-order weak approximations obtained from some implicit schemes. In [12], [13], we proposed constructing positive weak first- and second-order approximations by splitting the deterministic and stochastic parts of a stochastic differential equation (SDE) (re)written in the Stratonovich form. Independently, Ninomiya and Victoir [14] used a similar idea in the multidimensional case. Unfortunately, in the case of CIR equation, the corresponding split-step approximations are not well defined for large $\sigma$ (when ${\sigma }^2>4a$). Alfonsi [3] solves this problem by switching to another scheme in a neighborhood of zero.

In this paper, we propose to construct weak second-order approximations of the CIR process by splitting the CIR equation into the deterministic part and the stochastic part in the Itô form, instead of the Stratonovich form.