--- a/financepy/models/cir_montecarlo.py
+++ b/financepy/models/cir_montecarlo.py
@@ -102,21 +102,21 @@
     if t == 0.0:
         return 1.0
 
-    if sigma == 0.0:
-        integral_mean = b * t + (r0 - b) * (1.0 - np.exp(-a * t)) / a
-        return np.exp(-integral_mean)
-
-    sigma2 = sigma * sigma
-    h = np.sqrt(a * a + 2.0 * sigma2)
-    exp_ht = np.exp(h * t)
-
-    denom = 2.0 * h + (a + h) * (exp_ht - 1.0)
-
-    aa = (2.0 * h * np.exp((a + h) * t / 2.0) / denom) ** (2.0 * a * b / sigma2)
-
-    bb = 2.0 * (exp_ht - 1.0) / denom
-
-    return aa * np.exp(-r0 * bb)
+    # Factor the affine coefficients before evaluation. This removes the
+    # enormous 1/sigma**2 exponent and exp(h*t), and extends continuously
+    # to sigma=0 without an arbitrary volatility cutoff.
+    h = np.hypot(a, np.sqrt(2.0) * sigma)
+    u = -np.expm1(-h * t)
+    x = (sigma / h) * (sigma / (h + a)) * u
+
+    if x == 0.0:
+        log_ratio = 1.0
+    else:
+        log_ratio = -np.log1p(-x) / x
+
+    bb = (u / h) / (1.0 - x)
+    log_aa = (2.0 * a * b / (h + a)) * ((u / h) * log_ratio - t)
+    return np.exp(log_aa - r0 * bb)
 
 
 @njit(