CIR¶

The Cox–Ingersoll–Ross (CIR) model

quantflow.sp.cir.CIR `pydantic-model` ¶

Bases: IntensityProcess

The Cox-Ingersoll-Ross (CIR) model is a mean-reverting square-root diffusion process.

\[\begin{equation} dx_t = \kappa (\theta - x_t) dt + \sigma \sqrt{x_t}\, dw_t \end{equation}\]

where \(w_t\) is a standard Wiener process. The process stays strictly positive (Feller condition) when

\[\begin{equation} 2\kappa\theta \geq \sigma^2 \end{equation}\]

Fields:

rate (float)
kappa (float)
sigma (float)
theta (float)
sample_algo (SamplingAlgorithm)

sigma `pydantic-field` ¶

sigma = 1.0

Volatility \(\sigma\)

theta `pydantic-field` ¶

theta = 1.0

Mean rate \(\theta\)

sample_algo `pydantic-field` ¶

sample_algo = IMPLICIT

Sampling algorithm

sigma2 `property` ¶

sigma2

The square of the volatility parameter, used in various calculations.

feller_condition `property` ¶

feller_condition

Value of \(2\kappa\theta - \sigma^2\); positive means the Feller condition holds.

is_positive `property` ¶

is_positive

True if the Feller condition holds, guaranteeing the process stays strictly positive.

rate `pydantic-field` ¶

rate = 1.0

Instantaneous initial rate \(x_0\)

kappa `pydantic-field` ¶

kappa = 1.0

Mean reversion speed \(\kappa\)

sample ¶

sample(paths, time_horizon=1, time_steps=100)

Source code in quantflow/sp/cir.py

def sample(
    self, paths: int, time_horizon: float = 1, time_steps: int = 100
) -> Paths:
    draws = Paths.normal_draws(paths, time_horizon, time_steps)
    return self.sample_from_draws(draws)

sample_from_draws ¶

sample_from_draws(paths, *args)

Source code in quantflow/sp/cir.py

def sample_from_draws(self, paths: Paths, *args: Paths) -> Paths:
    match self.sample_algo:
        case SamplingAlgorithm.EULER:
            return self.sample_euler(paths)
        case SamplingAlgorithm.MILSTEIN:
            return self.sample_euler(paths, 0.25)
        case SamplingAlgorithm.IMPLICIT:
            return self.sample_implicit(paths)

sample_euler ¶

sample_euler(draws, milstein_coef=0.0)

Source code in quantflow/sp/cir.py

def sample_euler(self, draws: Paths, milstein_coef: float = 0.0) -> Paths:
    kappa = self.kappa
    theta = self.theta
    dt = draws.dt
    sdt = self.sigma * np.sqrt(dt)
    sdt2 = sdt * sdt
    paths = np.zeros(draws.data.shape)
    paths[0, :] = self.rate
    for t in range(draws.time_steps):
        w = sdt * draws.data[t, :]
        x = paths[t, :]
        xplus = np.clip(x, 0, None)
        dx = (
            kappa * (theta - xplus) * dt
            + np.sqrt(xplus) * w
            + milstein_coef * (w * w - sdt2)
        )
        paths[t + 1, :] = x + dx
    return Paths(t=draws.t, data=paths)

sample_implicit ¶

sample_implicit(draws)

Use an implicit scheme to preserve positivity of the process.

Source code in quantflow/sp/cir.py

def sample_implicit(self, draws: Paths) -> Paths:
    """Use an implicit scheme to preserve positivity of the process."""
    kappa = self.kappa
    theta = self.theta
    dt = draws.dt
    kdt2 = 2 * (1 + kappa * dt)
    kts = (kappa * theta - 0.5 * self.sigma2) * dt
    sdt = self.sigma * np.sqrt(dt)
    paths = np.zeros(draws.data.shape)
    paths[0, :] = self.rate
    for t in range(draws.time_steps):
        w = sdt * draws.data[t, :]
        x = paths[t, :]
        w2p = np.clip(w * w + 2 * (x + kts) * kdt2, 0, None)
        xs = (w + np.sqrt(w2p)) / kdt2
        paths[t + 1, :] = xs * xs
    return Paths(t=draws.t, data=paths)

characteristic_exponent ¶

characteristic_exponent(t, u)

Characteristic exponent of the CIR process.

\[\begin{equation} \begin{aligned} \phi(t, u) &= -\frac{2\kappa\theta}{\sigma^2} \!\left(\kappa t + \log\frac{2\kappa}{c}\right) - \frac{2\kappa\, iu}{c}\, x_0 \\ c &= iu\sigma^2 + (2\kappa - iu\sigma^2)\,e^{\kappa t} \end{aligned} \end{equation}\]

where \(x_0\) is the initial rate.

Source code in quantflow/sp/cir.py

def characteristic_exponent(self, t: Vector, u: Vector) -> Vector:
    r"""Characteristic exponent of the CIR process.

    \begin{equation}
    \begin{aligned}
        \phi(t, u) &= -\frac{2\kappa\theta}{\sigma^2}
        \!\left(\kappa t + \log\frac{2\kappa}{c}\right)
        - \frac{2\kappa\, iu}{c}\, x_0 \\
        c &= iu\sigma^2 + (2\kappa - iu\sigma^2)\,e^{\kappa t}
    \end{aligned}
    \end{equation}

    where $x_0$ is the initial rate.
    """
    iu = 1j * u
    sigma = self.sigma
    kappa = self.kappa
    kt = kappa * t
    ekt = np.exp(kt)
    sigma2 = sigma * sigma
    s2u = iu * sigma2
    c = s2u + (2 * kappa - s2u) * ekt
    b = 2 * kappa * iu / c
    a = 2 * kappa * self.theta * (kt + np.log(2 * kappa / c)) / sigma2
    return -a - b * self.rate

integrated_log_laplace ¶

integrated_log_laplace(t, u)

Log-Laplace transform of the time-integrated CIR process.

\[\begin{equation} \phi(t, u) = \log E\!\left[e^{-u \int_0^t x_s\, ds}\right] = \frac{2\kappa\theta}{\sigma^2} \log\!\left(\frac{2\gamma\, e^{(\gamma+\kappa)t/2}}{D}\right) - \frac{2u(e^{\gamma t}-1)}{D}\, x_0 \end{equation}\]

where \(\gamma = \sqrt{\kappa^2 + 2u\sigma^2}\), \(D = 2\gamma + (\gamma+\kappa)(e^{\gamma t}-1)\), and \(x_0\) is the initial rate.

Source code in quantflow/sp/cir.py

def integrated_log_laplace(self, t: Vector, u: Vector) -> Vector:
    r"""Log-Laplace transform of the time-integrated CIR process.

    \begin{equation}
        \phi(t, u) = \log E\!\left[e^{-u \int_0^t x_s\, ds}\right]
        = \frac{2\kappa\theta}{\sigma^2}
          \log\!\left(\frac{2\gamma\, e^{(\gamma+\kappa)t/2}}{D}\right)
        - \frac{2u(e^{\gamma t}-1)}{D}\, x_0
    \end{equation}

    where $\gamma = \sqrt{\kappa^2 + 2u\sigma^2}$,
    $D = 2\gamma + (\gamma+\kappa)(e^{\gamma t}-1)$, and $x_0$ is
    the initial rate.
    """
    kappa = self.kappa
    sigma2 = self.sigma2
    gamma = np.sqrt(kappa * kappa + 2 * u * sigma2)
    kts = 2 * kappa * self.theta / sigma2
    gamma_kappa = gamma + kappa
    egt = np.exp(gamma * t)
    d = 2 * gamma + gamma_kappa * (egt - 1)
    a = 2 * gamma * np.exp(0.5 * gamma_kappa * t) / d
    b = 2 * u * (egt - 1) / d
    return kts * np.log(a) - b * self.rate

domain_range ¶

domain_range()

Source code in quantflow/sp/cir.py

def domain_range(self) -> Bounds:
    return Bounds(0, np.inf)

analytical_mean ¶

analytical_mean(t)

Analytical mean of the CIR process at time \(t\).

\[\begin{equation} E[x_t] = x_0\, e^{-\kappa t} + \theta\bigl(1 - e^{-\kappa t}\bigr) \end{equation}\]

Source code in quantflow/sp/cir.py

def analytical_mean(self, t: FloatArrayLike) -> FloatArrayLike:
    r"""Analytical mean of the CIR process at time $t$.

    \begin{equation}
        E[x_t] = x_0\, e^{-\kappa t} + \theta\bigl(1 - e^{-\kappa t}\bigr)
    \end{equation}
    """
    ekt = self.ekt(t)
    return self.rate * ekt + self.theta * (1 - ekt)

analytical_variance ¶

analytical_variance(t)

Analytical variance of the CIR process at time \(t\).

\[\begin{equation} \mathrm{Var}(x_t) = \frac{\sigma^2(1 - e^{-\kappa t})}{\kappa} \left(x_0\, e^{-\kappa t} + \frac{\theta}{2}(1 - e^{-\kappa t})\right) \end{equation}\]

Source code in quantflow/sp/cir.py

def analytical_variance(self, t: FloatArrayLike) -> FloatArrayLike:
    r"""Analytical variance of the CIR process at time $t$.

    \begin{equation}
        \mathrm{Var}(x_t) = \frac{\sigma^2(1 - e^{-\kappa t})}{\kappa}
        \left(x_0\, e^{-\kappa t} + \frac{\theta}{2}(1 - e^{-\kappa t})\right)
    \end{equation}
    """
    kappa = self.kappa
    ekt = self.ekt(t)
    return (
        self.sigma2
        * (1 - ekt)
        * (self.rate * ekt + 0.5 * self.theta * (1 - ekt))
        / kappa
    )

analytical_pdf ¶

analytical_pdf(t, x)

The marginal pdf of the CIR process is the scaled non-central chi-squared.

\[\begin{equation} \begin{aligned} p(x_t = x \mid x_0) &= c\, e^{-(u+v)} \left(\frac{v}{u}\right)^{q/2} I_q\!\left(2\sqrt{uv}\right) \\ c &= \frac{2\kappa}{\sigma^2(1 - e^{-\kappa t})} \\ u &= c\, x_0\, e^{-\kappa t} \\ v &= c\, x \\ q &= \frac{2\kappa\theta}{\sigma^2} - 1 \end{aligned} \end{equation}\]

\(I_q\) is the modified Bessel function of the first kind of order \(q\).

Source code in quantflow/sp/cir.py

def analytical_pdf(self, t: FloatArrayLike, x: FloatArrayLike) -> FloatArrayLike:
    r"""The marginal pdf of the CIR process is the scaled non-central chi-squared.

    \begin{equation}
    \begin{aligned}
        p(x_t = x \mid x_0) &= c\, e^{-(u+v)}
            \left(\frac{v}{u}\right)^{q/2} I_q\!\left(2\sqrt{uv}\right) \\
        c &= \frac{2\kappa}{\sigma^2(1 - e^{-\kappa t})} \\
        u &= c\, x_0\, e^{-\kappa t} \\
        v &= c\, x \\
        q &= \frac{2\kappa\theta}{\sigma^2} - 1
    \end{aligned}
    \end{equation}

    $I_q$ is the modified Bessel function of the first kind of order $q$.
    """
    k = self.kappa
    s2 = self.sigma2
    ekt = self.ekt(t)
    c = 2 * k / (1 - ekt) / s2
    q = 2 * k * self.theta / s2 - 1
    u = c * ekt * self.rate
    v = c * x
    return (
        c
        * np.exp(-v - u)
        * np.power(v / u, 0.5 * q)
        * special.iv(q, 2 * np.sqrt(u * v))
    )

characteristic ¶

characteristic(t, u)

Characteristic function at time t for a given input parameter u

The characteristic function represents the Fourier transform of the probability density function

\[\begin{equation} \Phi = {\mathbb E} \left[e^{i u x_t}\right] = e^{-\phi(t, u)} \end{equation}\]

where \(\phi\) is the characteristic exponent, which can be more easily computed for many processes.

PARAMETER	DESCRIPTION
`t`	Time horizon TYPE: `FloatArrayLike`
`u`	Characteristic function input parameter TYPE: `Vector`

Source code in quantflow/sp/base.py

def characteristic(
    self,
    t: Annotated[FloatArrayLike, Doc("Time horizon")],
    u: Annotated[Vector, Doc("Characteristic function input parameter")],
) -> Vector:
    r"""Characteristic function at time `t` for a given input parameter `u`

    The characteristic function represents the Fourier transform of the
    probability density function

    \begin{equation}
        \Phi = {\mathbb E} \left[e^{i u x_t}\right] = e^{-\phi(t, u)}
    \end{equation}

    where $\phi$ is the characteristic exponent, which can be more easily
    computed for many processes.
    """
    return np.exp(-self.characteristic_exponent(t, u))

convexity_correction ¶

convexity_correction(t)

Convexity correction for the process

Source code in quantflow/sp/base.py

def convexity_correction(self, t: FloatArrayLike) -> Vector:
    """Convexity correction for the process"""
    return -self.characteristic_exponent(t, complex(0, -1)).real

analytical_std ¶

analytical_std(t)

Analytical standard deviation of the process at time t

This has a closed form solution if the process has an analytical variance

Source code in quantflow/sp/base.py

def analytical_std(self, t: FloatArrayLike) -> FloatArrayLike:
    """Analytical standard deviation of the process at time `t`

    This has a closed form solution if the process has an analytical variance
    """
    return np.sqrt(self.analytical_variance(t))

analytical_cdf ¶

analytical_cdf(t, x)

Analytical cdf of the process at time t

Implement if available

Source code in quantflow/sp/base.py

def analytical_cdf(self, t: FloatArrayLike, x: FloatArrayLike) -> FloatArrayLike:
    """Analytical cdf of the process at time `t`

    Implement if available
    """
    raise NotImplementedError

marginal ¶

marginal(t)

Marginal distribution of the process at time t

Source code in quantflow/sp/base.py

def marginal(self, t: FloatArrayLike) -> StochasticProcess1DMarginal:
    """Marginal distribution of the process at time `t`"""
    return StochasticProcess1DMarginal(process=self, t=t)

frequency_range ¶

frequency_range(std, max_frequency=None)

Maximum frequency when calculating characteristic functions

Source code in quantflow/sp/base.py

def frequency_range(self, std: float, max_frequency: float | None = None) -> Bounds:
    """Maximum frequency when calculating characteristic functions"""
    if max_frequency is None:
        max_frequency = np.sqrt(40 / std / std)
    return Bounds(0, max_frequency)

support ¶

support(mean, std, points)

Support of the process at time t

Source code in quantflow/sp/base.py

def support(self, mean: float, std: float, points: int) -> FloatArray:
    """Support of the process at time `t`"""
    bounds = self.domain_range()
    start = float(sigfig(bound_from_any(bounds.lb, mean - std)))
    end = float(sigfig(bound_from_any(bounds.ub, mean + std)))
    return np.linspace(start, end, points + 1)

ekt ¶

ekt(t)

Source code in quantflow/sp/base.py

def ekt(self, t: FloatArrayLike) -> FloatArrayLike:
    return np.exp(-self.kappa * t)

CIR¶

quantflow.sp.cir.CIR pydantic-model ¶

sigma pydantic-field ¶

theta pydantic-field ¶

sample_algo pydantic-field ¶

sigma2 property ¶

feller_condition property ¶

is_positive property ¶

rate pydantic-field ¶

kappa pydantic-field ¶

sample ¶

sample_from_draws ¶

sample_euler ¶

sample_implicit ¶

characteristic_exponent ¶

integrated_log_laplace ¶

domain_range ¶

analytical_mean ¶

analytical_variance ¶

analytical_pdf ¶

characteristic ¶

convexity_correction ¶

analytical_std ¶

analytical_cdf ¶

marginal ¶

frequency_range ¶

support ¶

ekt ¶

quantflow.sp.cir.CIR `pydantic-model` ¶

sigma `pydantic-field` ¶

theta `pydantic-field` ¶

sample_algo `pydantic-field` ¶

sigma2 `property` ¶

feller_condition `property` ¶

is_positive `property` ¶

rate `pydantic-field` ¶

kappa `pydantic-field` ¶