Ornstein-Uhlenbeck Density from the Fokker-Planck Equation

From Mean Reversion to a Time-Dependent Density

An Ornstein–Uhlenbeck model describes a fluctuation that is pulled back toward a reference level while being continuously perturbed by random noise. In pharmacokinetics, the state $X_t$ can represent a concentration fluctuation around a typical concentration $\mu$. A fluctuation is not itself a total concentration: it may be negative without implying that the corresponding total concentration is negative.

The model has linear mean-reverting drift and constant diffusion. The forward equation describes how the probability density changes over time. Here the goal is to solve that equation and determine the full distribution, not only its first two moments.

Definition: Let $X_t$ satisfy the Ornstein–Uhlenbeck stochastic differential equation $$ dX_t=-\lambda(X_t-\mu)\,dt+\sigma\,dW_t, $$ where $\lambda>0$ is the mean-reversion rate, $\mu$ is the reference level, $\sigma>0$ is the diffusion coefficient, and $W_t$ is a Wiener process. Assume $X_0$ is independent of $W_t$. The forward equation for a density $p(x,t)$ is $$ \frac{\partial p}{\partial t} =\lambda\frac{\partial}{\partial x}\bigl((x-\mu)p\bigr) +\frac{\sigma^2}{2}\frac{\partial^2p}{\partial x^2}. $$

The coefficient of the second derivative is $\sigma^2/2$, not $\sigma/2$: the forward equation uses the diffusion variance per unit time. In the derivation below, assume the density and its derivatives decay sufficiently rapidly as $|x|\to\infty$ for the indicated integrations by parts to have vanishing boundary terms. The result also applies to a deterministic initial state, interpreted as a point-mass initial distribution; for $\sigma>0$, the density is smooth for every $t>0$.

Deriving the Gaussian Law from the Forward Equation

Center the state at its reference level by setting $Y_t=X_t-\mu$. Define its characteristic function by $\psi(u,t)=\int_{\mathbb R}e^{iuy}p_Y(y,t)\,dy$. The centered forward equation is $\partial_t p_Y=\lambda\partial_y(yp_Y)+(\sigma^2/2)\partial_{yy}p_Y$. Taking its Fourier transform gives $\partial_t\psi=-\lambda u\,\partial_u\psi-(\sigma^2u^2/2)\psi$. Indeed, integration by parts gives $\int e^{iuy}\partial_y(yp_Y)\,dy=-u\,\partial_u\psi$ and $\int e^{iuy}\partial_{yy}p_Y\,dy=-u^2\psi$.

This first-order equation can be solved along characteristics. Fix a terminal frequency $u$ at time $t$. The frequency along the characteristic at time $s$ is $u(s)=u e^{-\lambda(t-s)}$. Along that path, $\frac{d}{ds}\psi(u(s),s)=-(\sigma^2/2)u(s)^2\psi(u(s),s)$. Integrating this scalar equation from $s=0$ to $s=t$ yields $\psi(u,t)=\psi(u e^{-\lambda t},0)\exp[-\frac{\sigma^2u^2}{2}\int_0^t e^{-2\lambda(t-s)}ds]$. The integral is $(1-e^{-2\lambda t})/(2\lambda)$.

Theorem: Suppose $X_0\sim N(m_0,v_0)$, where $v_0\geq 0$, and $X_0$ is independent of $W_t$. For $t>0$, the solution of the Ornstein–Uhlenbeck forward equation has Gaussian distribution with $$ m(t)=E[X_t]=\mu+(m_0-\mu)e^{-\lambda t}, \qquad V(t)=\operatorname{Var}(X_t) =v_0e^{-2\lambda t}+\frac{\sigma^2}{2\lambda}(1-e^{-2\lambda t}). $$ Thus $X_t\sim N(m(t),V(t))$. When $v_0=0$, the initial distribution is a point mass at $m_0$, and the stated density applies for $t>0$.

To establish the theorem, the centered initial characteristic function is $\psi(u,0)=\exp(iu(m_0-\mu)-v_0u^2/2)$. Substituting $u e^{-\lambda t}$ for its frequency and evaluating the characteristic-path integral gives $$ \psi(u,t)= \exp\left(iu(m_0-\mu)e^{-\lambda t} -\frac{u^2}{2}\left[v_0e^{-2\lambda t} +\frac{\sigma^2}{2\lambda}(1-e^{-2\lambda t})\right]\right). $$ This is the characteristic function of a Gaussian variable with the displayed centered mean and variance. Shifting back from $Y_t$ to $X_t$ adds $\mu$ to the mean and leaves the variance unchanged. Characteristic functions uniquely determine probability distributions, so this identifies the law and the corresponding density. In particular, for $V(t)>0$, $$ p(x,t)=\frac{1}{\sqrt{2\pi V(t)}} \exp\left(-\frac{(x-m(t))^2}{2V(t)}\right). $$ $\square$

The derivation solves the forward equation by transforming it into an equation for the characteristic function. The Gaussian conclusion follows from the explicit solution of that equation, not merely from matching the first two moments. This distinction matters: knowing a mean and variance alone does not determine a general probability distribution.

Worked Densities and Checks

Worked Example: A Noisy Fluctuation After One Hour

Suppose $\lambda=0.5\ \mathrm{h}^{-1}$, $\mu=0.30\ \mathrm{mg/L}$, $X_0\sim N(0.50,0.040)$, and $\sigma^2=0.090\ (\mathrm{mg/L})^2/\mathrm{h}$. Find the distribution after one hour.

The mean is $m(1)=0.30+(0.50-0.30)e^{-0.5} =0.30+0.20(0.6065306597) =0.4213061319\ \mathrm{mg/L}$. For the variance, $2\lambda=1\ \mathrm{h}^{-1}$, so $$ V(1)=0.040e^{-1} +\frac{0.090}{1}(1-e^{-1}) =0.0147151776+0.0568908503 =0.0716060279\ (\mathrm{mg/L})^2. $$

The correctly rounded variance is $0.0716060\ (\mathrm{mg/L})^2$ to seven decimal places. Therefore $X_1\sim N(0.421306,\ 0.0716060)$, with the second parameter denoting variance. The density is centered above $\mu$, as expected because the starting mean $0.50$ exceeds the reference level and has not yet fully returned to it.

Worked Example: Transition Density from a Fixed Starting Value

Let $X_0=0.10\ \mathrm{mg/L}$, $\mu=0.50\ \mathrm{mg/L}$, $\lambda=1\ \mathrm{h}^{-1}$, and $\sigma^2=0.080\ (\mathrm{mg/L})^2/\mathrm{h}$. Find the transition density at $t=\ln 2$ hours.

Here $e^{-\lambda t}=e^{-\ln 2}=1/2$. Since the starting value is fixed, $v_0=0$, and the mean is $m(t)=0.50+(0.10-0.50)(1/2)=0.30\ \mathrm{mg/L}$. The variance is $$ V(t)=\frac{0.080}{2(1)}\left(1-e^{-2\ln 2}\right) =0.040\left(1-\frac14\right) =0.030\ (\mathrm{mg/L})^2. $$

Thus the transition density is $N(0.30,0.030)$. Although the initial state had no uncertainty, the positive diffusion coefficient creates a nonzero variance at every positive time.

Worked Example: Checking a Stationary Gaussian Distribution

Suppose $X_0\sim N(0.60,0.040)$, $\mu=0.60\ \mathrm{mg/L}$, $\lambda=0.25\ \mathrm{h}^{-1}$, and $\sigma^2=0.020\ (\mathrm{mg/L})^2/\mathrm{h}$. Check the distribution at $t=3$ hours.

Since $m_0=\mu$, the mean remains $m(3)=0.60\ \mathrm{mg/L}$. Also, $\sigma^2/(2\lambda)=0.020/0.50=0.040\ (\mathrm{mg/L})^2$ and $e^{-2\lambda t}=e^{-1.5}$. Therefore $$ V(3)=0.040e^{-1.5}+0.040(1-e^{-1.5}) =0.040\ (\mathrm{mg/L})^2. $$

The law remains $N(0.60,0.040)$. More generally, the Gaussian law with mean $\mu$ and variance $\sigma^2/(2\lambda)$ is invariant: inserting those values as $m_0$ and $v_0$ leaves both the mean and variance unchanged for every $t\geq0$.

Assumptions, Diagnostics, and Alternatives

The formula requires linear mean reversion, constant diffusion, and an initial Gaussian law if the marginal distribution itself is to be Gaussian at every time. With a non-Gaussian initial distribution, the conditional distribution given $X_0=x_0$ is still Gaussian, but the unconditional law is generally a mixture of such Gaussians and need not be Gaussian. The formula also presumes an unconstrained real-valued state; it is a model for a fluctuation, not by itself a guarantee that a total concentration remains nonnegative.

Several checks can reveal errors. As $t$ approaches zero, $m(t)$ must approach $m_0$ and $V(t)$ must approach $v_0$. As $t$ grows, the mean approaches $\mu$, while the variance approaches $\sigma^2/(2\lambda)$. The variance must be nonnegative. Its noise contribution must increase toward that limit, while the contribution from the initial variance decays at rate $2\lambda$, twice the exponent rate for the mean.

A frequent mistake is to use $e^{-\lambda t}$ instead of $e^{-2\lambda t}$ in the variance formula, or to use $\sigma$ where $\sigma^2$ is required. Another is to round separate variance components and then report their rounded sum as though it were the correctly rounded total. In the first example, the unrounded total is $0.0716060279\ldots$; rounding the two displayed components first would instead produce a slightly different final digit. Keep precision through the addition, then round once.

The same law can be derived by solving the linear SDE with an integrating factor and evaluating the mean and variance of its stochastic integral. That route is often efficient for moments and makes the effect of the initial state transparent. The forward-equation route instead directly solves for the density and generalizes conceptually to cases where a density equation is the main object of interest. Neither route removes the model assumptions: if drift or diffusion changes with state, this Gaussian solution does not automatically apply.

Key takeaway: For a Gaussian initial state independent of the Wiener process, an Ornstein–Uhlenbeck fluctuation remains Gaussian. Its mean relaxes toward $\mu$ at rate $\lambda$, and its variance relaxes toward $\sigma^2/(2\lambda)$ at rate $2\lambda$.

Check Your Understanding

Answer each question, then compare with the solutions below.

If $m_0=\mu$, what is $m(t)$, and why?
For a deterministic start, write the variance at time $t>0$.
With $\lambda=0.5\ \mathrm{h}^{-1}$ and $\sigma^2=0.10\ (\mathrm{mg/L})^2/\mathrm{h}$, what is the limiting variance?
Does a Gaussian transition density conditional on $X_0$ guarantee a Gaussian marginal density for a non-Gaussian $X_0$? Explain.

Solutions. (1) $m(t)=\mu+(\mu-\mu)e^{-\lambda t}=\mu$. (2) Set $v_0=0$, giving $V(t)=\frac{\sigma^2}{2\lambda}(1-e^{-2\lambda t})$. (3) The limit is $\sigma^2/(2\lambda)=0.10/1.0=0.10\ (\mathrm{mg/L})^2$. (4) No. Averaging the conditional Gaussian laws over a non-Gaussian initial distribution generally gives a Gaussian mixture, which need not be Gaussian.