A square-root diffusion process $X$—also known a Cox-Ingersoll-Ross (CIR) process—is that which is defined by the following stochastic differential equation

$$ {\rm d}X_t = a\big( b - X_t \big) {\rm d}t + \sigma\sqrt{X_t} {\rm d}W_t, $$

where $W$ is a standard Brownian motion on $\bbR$ and $a, b, \sigma \in \bbR$ satsify $2ab \geq \sigma^2$. Such a process is nice to play with analytically, as it has very tractable properties; most notably, it is time-homogeneous, Markov, affine, and ergodic. In this post, I show how the affine structure of a CIR process $X$ allows us to easily calculate its long-time distribution $\mu_X$. We then demonstrate numerically how samples $\rmP(X_T \in \cdot | X_0 = x)$ distribute like $\mu_X$ for large times $T$.

Long-time distribution

Definition. The long-time distribution of a stochastic process $X$ on $\bbR$ should be a measure $\mu_X$ satisfying the following weak limit.

$$ \lim_{t\rightarrow\infty} \rmP(X_t \in \cdot) = \mu_X $$

A classic result of Lévy (see Kallenberg[undefined] or Wikipedia[undefined]) dictates that such a weak convergence is determined by that of pointwise convergence of the characteristic functions.

$$ \lim_{t\rightarrow\infty} \rmE\big( e^{iu X_t} \big) = \int e^{iuy} \mu_X({\rm d} y), \quad u \in \bbR $$

Because the drift $x \mapsto a(b - x)$ and the square-diffusion $x \mapsto \sigma^2 x$ of our CIR process $X$ are affine functions, we are able to conclude that $X$ is an affine process (see this post on affine processes). Thus, solving the expectation $\rm E e^{iu X_t}$ for each $u \in \bbR$ and $t \geq 0$ amounts to evaluating

$$ \rmE \big( e^{iu X_t} | X_0 = x \big) = e^{\psi_0(t, u) + \psi_1(t, u)x} $$

for $\psi_0, \psi_1$ solving the following differential equation.

A simple separation of variables allows us to solve the above ordinary differential equation; from here, the characteristic function is as follows.

$$ \rmE\big( e^{iu X_t} | X_0 = x \big) = \exp\Bigg(\frac{2ab}{\sigma^2}\log\frac{2ae^{at}}{i\sigma^2u + (2a-i\sigma^2 u)e^{at}} + \frac{i2au}{i\sigma^2u + (2a - i\sigma^2u)e^{at}} x \Bigg) $$

It is now clear that our pointwise limit is as follows.

$$ \lim_{t\rightarrow\infty} \rmE\big( e^{iu X_t} | X_0 = x \big) = \exp\Bigg(\frac{2ab}{\sigma^2}\log\Big(\frac{2a}{2a-i\sigma^2 u}\Big)\Bigg) = \Big( 1 - \frac{\sigma^2}{2a}iu\Big)^{-2ab/\sigma^2} $$

The keen eye will recognize that the final expression above is the characteristic function of a Gamma distribution[undefined] of shape $2ab/\sigma^2$ and rate $2a/\sigma^2$.

$$ \mu_X = \operatorname{Gamma}\Big(\text{shape} = \frac{2ab}{\sigma^2}, \text{rate} = \frac{2a}{\sigma^2} \Big) $$

Its density and moment generating functions have closed forms, up to evaluation of the Gamma function.

Numerics

From the perspective of Monte Carlo, the fact that $X$ has a long-time distribution $\mu_X$ should mean that samples of a marginal $X_T$ for a large time $T$ should distribute similar to $\mu_X$. We intend to show this by taking such samples and comparing the empirical moment generating function and density to those of $\mu_X$ above.

Below, we write small Julia[undefined] script which invokes a simple Euler-Maruyama[undefined] scheme to simulate $100$ (biased) samples of marginals

import Random: randn
import JSON: json

# CIR process with parameters (a, b, σ, X_0) = (1.2, 3.1, 0.5, 2.2)
a = 1.2
b = 3.1
sigma = 0.5
x0 = 2.2

# generate samples
T = 10.0
times = LinRange(0.0, T, 2501)
count = 100
trajectories = zeros(length(times), count)
normals = randn(length(times)-1, count)
for i in 1:count
  trajectories[1,i] = x0
  for k in 2:length(times)
    dt = times[k] - times[k-1]
    dw = normals[k-1, i] * sqrt(dt)
    x = trajectories[k-1,i]
    trajectories[k,i] = x + a * (b - x) * dt + sigma * sqrt(abs(x)) * dw
  end
end

# export to json
write(
  "data.json",
  json(Dict(
    :a => a, 
    :b => b, 
    :sigma => sigma, 
    :samples => map(
      i -> map(
        ((x, t), ) -> Dict(:xt => x, :t => t, :sample => i), 
        zip(trajectories[:,i], times)
      ),
      1:count
    )
  ))
)

With these samples, we are able to visualize the trajectories and statistics of the empirical distribution. In the visualization below, the $100$ blue paths demonstrate our sample trajectories $t \mapsto X_t$. The green dots represent the empirical distribution of $X_t$ for the time $t$ over which the mouse currently hovers. For such an empirical distribution, we may compare the moment generating function and a kernel density estimation of the density to those of the real long-term distribution $\mu_X$.

Takeaway

From the perspective of a mathematician, the closed form of the characteristic function $$u \mapsto \rmE(e^{iuX_t} | X_0 = x),$$ or similarly moment generating function $$u \mapsto \rmE(e^{uX_t} | X_0 = x),$$ helps make analytical statements regarding the process $X$. From the perspective of someone interested in computer simulation, costly samples of $X_T$ could potentially be approximated by efficient samples of $\mu_X$.

References

Related posts