Why Does CIR Go Negative in Simulation? (And How to Fix It)

You implemented the CIR model, ran a simulation, and your rates went negative — or worse, your code threw a wall of NaNs. If you’re pricing with a square-root process (CIR for rates, Heston for variance) and hitting this, you’ve run into one of the most common and most misunderstood bugs in quant implementation.

The short version: it’s not really a bug in your code, it’s a mismatch between the continuous model and the discrete scheme you’re simulating it with. Here’s what’s happening and how to fix it properly.

The model, and why it should stay positive

The CIR process is:

$dr_t = \kappa(\theta - r_t)\,dt + \sigma\sqrt{r_t}\,dW_t$

The whole point of the square root is to keep the process positive. As $r_t$ approaches zero, the $\sqrt{r_t}$ term shrinks the volatility to zero, so the random shocks fade out and the mean-reverting drift $\kappa(\theta - r_t)$ — which is positive when $r_t < \theta$ — pulls the process back up. In the continuous model, the diffusion switches off exactly when it would otherwise push you below zero. That’s the theory. It’s why CIR is used for quantities that can’t be negative: short rates, default intensities, stochastic variance.

The Feller condition: when zero is truly unreachable

The continuous process stays strictly positive only if the parameters satisfy the Feller condition:

$2\kappa\theta \geq \sigma^2$

Read it as a tug of war. The left side is the strength of the upward pull near zero (mean reversion times the level it reverts to). The right side is the strength of the volatility trying to push you through. If mean reversion wins, zero is unattainable and the process never touches it. If volatility wins, the process can reach zero — it bounces off rather than going negative, but it does get there.

kappa, theta, sigma = 0.5, 0.04, 0.30

print(f"2*kappa*theta = {2*kappa*theta:.4f}")
print(f"sigma^2       = {sigma**2:.4f}")
print(f"Feller holds? {2*kappa*theta >= sigma**2}")

2*kappa*theta = 0.0400
sigma^2       = 0.0900
Feller holds? False

With these parameters — a high volatility of 30% — Feller fails badly. Zero is attainable. And this is the regime where naive simulation falls apart, because calibrated Heston and CIR parameters routinely violate Feller. It’s not an edge case; it’s the normal case in practice.

Why the simulation goes negative when the model doesn’t

Here’s the crux. The continuous process never goes negative. Your simulation does. Why?

Because you’re not simulating the continuous process — you’re simulating a discrete approximation of it, and the approximation takes finite steps. The standard Euler scheme updates like this:

$r_{n+1} = r_n + \kappa(\theta - r_n)\Delta t + \sigma\sqrt{r_n}\sqrt{\Delta t}\,Z_n$

When $r_n$ is small and $Z_n$ is a large negative draw, that last term can be bigger than $r_n$ itself. The step overshoots zero and lands on a negative number — something the continuous process, which adjusts its volatility continuously, would never do. The discrete scheme only checks the volatility at the start of the step, so it can take too big a jump before the shrinking-volatility effect kicks in.

And then it gets worse on the very next step:

import numpy as np
np.random.seed(0)

kappa, theta, sigma, r0, T = 0.5, 0.04, 0.30, 0.04, 5.0
n_paths, n_steps = 20_000, 500
dt = T/n_steps

r = np.full(n_paths, r0)
for i in range(n_steps):
    Z = np.random.standard_normal(n_paths)
    r = r + kappa*(theta-r)*dt + sigma*np.sqrt(r)*np.sqrt(dt)*Z
    if np.isnan(r).any():
        print(f"NaN at step {i}: sqrt of a negative rate")
        break

NaN at step 217: sqrt of a negative rate

Once a rate goes negative, the next step needs $\sqrt{r_n}$ of a negative number, which is NaN. That NaN then contaminates everything downstream. The wall of NaNs you saw isn’t random — it’s one negative rate poisoning the rest of the simulation.

The fix: full truncation

The standard, battle-tested fix is the full truncation scheme. The idea is simple: wherever the process appears inside the square root and the drift, replace $r_n$ with $\max(r_n, 0)$ . You let the state variable go slightly negative, but you never feed a negative number into the square root.

r = np.full(n_paths, r0)
for i in range(n_steps):
    Z = np.random.standard_normal(n_paths)
    r_plus = np.maximum(r, 0.0)                      # truncate for the dynamics
    r = r + kappa*(theta - r_plus)*dt \
          + sigma*np.sqrt(r_plus)*np.sqrt(dt)*Z

Run this with the same Feller-violating parameters and it completes cleanly, with the long-run mean landing where theory says it should:

Final mean rate: 0.0401
Long-run mean theta: 0.04

No NaNs. The simulation runs, and the mean matches the theoretical $\theta = 0.04$ . Full truncation is the scheme most production systems use because it’s simple, robust, and has the best convergence behaviour among the naive fixes.

Fixes that seem reasonable but aren’t

People reach for these first, and they’re worse:

Taking the absolute value, $\sqrt{|r_n|}$ . It avoids the NaN, but it reflects the process off zero with the wrong dynamics and biases your results. Don’t.
Clamping the rate to a small floor like $10^{-8}$ . Stops the NaN but distorts the behaviour near zero — precisely the region you cared about modelling correctly. It also injects a systematic bias.
Just using more time steps. Helps, because smaller steps overshoot less often, but it doesn’t eliminate the problem and it’s expensive. You’ll still hit negative values occasionally, and “occasionally” is enough to produce a NaN that ruins a path.

Full truncation is preferred over all of these because it’s unbiased in the right way and cheap.

When you need it exact: the QE scheme

For cases where even full truncation’s small bias matters — long-dated exotics, or Heston calibration where accuracy near zero variance is critical — the gold standard is Andersen’s Quadratic-Exponential (QE) scheme. It samples from an approximation of the true non-central chi-squared distribution that CIR actually follows, rather than from a Gaussian step. It’s more work to implement, and it’s the right tool when full truncation isn’t accurate enough.

The CIR transition density is genuinely non-central chi-squared — you can even sample it exactly — but the exact method is slow, which is why practical schemes like full truncation and QE exist as the workhorses.

The takeaway

Your CIR simulation goes negative because the discrete Euler scheme takes finite steps that overshoot zero, even though the continuous process never would — and the problem is worst exactly when your calibrated parameters violate the Feller condition, which is most of the time. The fix is full truncation: apply $\max(r, 0)$ inside the square root and drift. For high-accuracy work, move to the QE scheme.

The broader lesson generalises to every SDE you simulate: the discretisation is not the model. Euler-Maruyama is convenient but naive, and knowing where it breaks — negative rates in CIR, the same issue in Heston’s variance — is a large part of what makes simulation code production-grade rather than merely plausible.

If you want to see this done properly — discretisation schemes, where they fail, and how to build simulation you can actually trust — it’s part of what we cover in Risk-Neutral Valuation, from the stochastic calculus up.

Get the Quant Reference Card — four pages covering stochastic calculus, measure change, the models you’ll actually meet, discretisation, rates, credit and XVA, with the practitioner notes that don’t make it into textbooks. Free.