The Square Root That Keeps Rain Positive: CIR Processes in Weather Modelling

From Finance to Meteorology

The Cox-Ingersoll-Ross process was introduced in 1985 as a model for short-term interest rates. Its defining feature — that it never becomes negative — was useful for rates, which were assumed to be positive. Three decades later, meteorologists and climate scientists discovered that the same property makes CIR a natural model for quantities like rainfall intensity, wind speed, and temperature variance: physical quantities that obey the same constraint for completely different reasons.

The migration of a stochastic process from fixed-income modelling to atmospheric science is not a coincidence. It reflects a structural property of the CIR SDE that happens to match the physics of non-negative meteorological variables.

The Equation and Why the Square Root Matters

The CIR process $X_t$ satisfies:

where $\kappa > 0$ is the speed of mean reversion, $\theta > 0$ is the long-run mean, and $\sigma > 0$ is the volatility parameter. The critical difference from the Ornstein-Uhlenbeck process is the $\sqrt{X_t}$ factor in the diffusion term.

This has two consequences. First, the variance of increments scales with the level: when $X_t$ is large, fluctuations are large; when $X_t$ is near zero, the noise term $\sigma\sqrt{X_t}$ vanishes. The process becomes quieter as it approaches zero, which is exactly what physically bounded processes do. Rainfall intensity does not fluctuate wildly when there is almost no rain. Second — and this is the key — the noise vanishes precisely at the boundary, which means the drift $\kappa(\theta – X_t) = \kappa\theta > 0$ at $X_t = 0$ pushes the process away before it can cross. Under the right condition, $X_t$ stays strictly positive forever.

The Feller Condition

The boundary behaviour of the CIR process is governed by a single inequality known as the **Feller condition**:

When this holds, the boundary at zero is inaccessible: the process is repelled from zero and $X_t > 0$ almost surely for all $t > 0$, regardless of initial condition. When $2\kappa\theta < \sigma^2$, zero is accessible and the process can touch — but not cross — the boundary, reflecting immediately.

For rainfall modelling, the Feller condition gives a direct constraint on calibrated parameters: if $\sigma^2 > 2\kappa\theta$, the fitted model predicts that dry spells (zero rainfall) occur with positive probability even within the continuous-time formulation, which may or may not be desirable depending on the application. For wind speed, where sustained calm is possible, a reflecting boundary may be the right choice.

The Stationary Distribution

Unlike the Ornstein-Uhlenbeck process, whose stationary distribution is Gaussian, the CIR process has a **Gamma stationary distribution**:

This is not a modelling assumption — it follows directly from the Fokker-Planck equation associated with (1). The Gamma distribution is right-skewed and supported on $(0, \infty)$, consistent with observed rainfall distributions, which are heavy-tailed and strictly positive. Gaussian models systematically underestimate the probability of extreme rainfall events. The CIR stationary distribution does not.

The shape parameter $2\kappa\theta/\sigma^2$ controls the skewness: smaller values (close to 1) produce distributions with a pronounced right tail, appropriate for arid climates where extreme rainfall events dominate the distribution. Larger values push the distribution toward symmetric, appropriate for temperate climates with more regular precipitation.

Simulating Rainfall Intensity

Because the CIR process has a known transition distribution (a scaled non-central chi-squared), it can be simulated exactly without discretisation error — analogous to the exact simulation of the OU process.

Input: kappa, theta, sigma, X_0, T, N
Set dt = T / N
Set c = sigma^2 * (1 - exp(-kappa*dt)) / (2*kappa)
Set nu = 4*kappa*theta / sigma^2          (degrees of freedom)

For i = 1 to N:
    lam = X[i-1] * exp(-kappa*dt) / c     (non-centrality parameter)
    X[i] = c * ChiSquared(nu, lam)         (non-central chi-squared draw)

Output: path X[0], X[1], ..., X[N]

The draw from the non-central chi-squared can be computed as a Poisson mixture of central chi-squared variates. This is the standard exact method used in Monte Carlo pricing of CIR-based bond models, and it carries over directly to meteorological simulation.

A Practical Application: Daily Rainfall Intensity

Consider modelling the daily rainfall intensity $R_t$ (in mm/day) at a given station over a season. The CIR model is:

$$dR_t = \kappa(\bar{R} – R_t)\,dt + \sigma\sqrt{R_t}\,dW_t$$

where $\bar{R}$ is the seasonal mean daily rainfall. Parameters $\kappa$, $\sigma$ are estimated from historical data by matching the sample mean, variance, and autocorrelation of the observed series. The fitted Gamma stationary distribution (3) then provides a parametric model for the rainfall distribution that can be used for drought probability estimation, flood frequency analysis, and crop yield modelling.

The CIR model does not attempt to capture every feature of the rainfall process — it ignores the wet/dry day structure, seasonal non-stationarity, and spatial correlation. But within its scope, it gives a physically consistent, analytically tractable, and easily calibrated model. That combination is rare.

Takeaway

The CIR process crossed into meteorology because it encodes, through the $\sqrt{X_t}$ diffusion term, a physical principle that the Ornstein-Uhlenbeck process ignores: quantities that cannot be negative become less volatile as they approach zero. That single structural feature — captured by one square root — changes the stationary distribution from Gaussian to Gamma, the boundary behaviour from accessible to inaccessible, and the fit to observed rainfall data from poor to good.

Interested in applying these ideas to your work? Get in touch.