Common Errors in Linear SDE Moment Calculations

Why Linear SDE Moment Calculations Go Wrong

Linear drift makes the mean of a linear SDE relatively direct to calculate. It does not make the second moments behave like the moments of a deterministic linear system. The crucial difference is that Brownian increments have nonzero quadratic variation: even when the mean follows a simple linear equation, diffusion contributes to the evolution of second moments.

Three errors are especially consequential in pharmacokinetic calculations: dropping that diffusion contribution, reporting a raw second moment as if it were a covariance, and replacing state-dependent diffusion with a constant without changing the model. The equations below provide checks that expose all three.

Moment Equations to Check Against

Consider a \(d\)-dimensional state \(X(t)\) driven by \(q\) independent standard Brownian motions \(W_1,\ldots,W_q\):

dX(t)=(A X(t)+a)\,dt+\sum_{j=1}^{q}(G_jX(t)+g_j)\,dW_j(t).

Here \(A\) and \(G_j\) are \(d\times d\) matrices, and \(a\) and \(g_j\) are \(d\)-vectors. The drift and diffusion coefficients are affine functions of the state. Define the mean \(m(t)=E[X(t)]\), the raw second-moment matrix \(M(t)=E[X(t)X(t)^T]\), and the covariance matrix \(P(t)=E[(X(t)-m(t))(X(t)-m(t))^T]\). Thus \(P=M-mm^T\); in general, \(M\) is not the covariance.

Assumptions: The coefficients are constant and finite, the Brownian motions are independent, and \(X(0)\) is independent of them with \(E[|X(0)|^4]<\infty\). The affine coefficients are globally Lipschitz and have at most linear growth, so the solution exists without explosion and has finite fourth moments on every finite time interval. These conditions justify taking expectations in the Itô product calculation below.

The mean equation is \(m'=Am+a\). It is useful, but it is not the second-moment equation. For example, it says nothing by itself about the noise-driven increase in variance.

Theorem: Under the stated assumptions, the raw second moment and covariance satisfy

\begin{aligned} M'={}&AM+MA^T+a m^T+m a^T\\ &+\sum_{j=1}^{q}\left(G_jMG_j^T+G_jm g_j^T+g_jm^TG_j^T+g_jg_j^T\right),\\[4pt] P'={}&AP+PA^T+\sum_{j=1}^{q}\left[G_jPG_j^T+(G_jm+g_j)(G_jm+g_j)^T\right]. \end{aligned}

The sum in either equation is the diffusion contribution. It is present even if the drift is deterministic and linear. The covariance equation also makes clear that, when the diffusion depends on the state, its contribution depends on both the covariance and the mean.

Derivation: Keep the Quadratic Variation

Apply the Itô product rule to \(XX^T\). Write \(b_j(X)=G_jX+g_j\). Since the Brownian drivers are independent, their cross-variation terms vanish, and

d(XX^T)=dX\,X^T+X\,dX^T+\sum_{j=1}^{q}b_j(X)b_j(X)^T\,dt.

The last term is the quadratic-variation contribution that is lost if one treats the product rule as ordinary differentiation. Taking expectations gives the raw-moment equation: the drift terms become \(AM+MA^T+a m^T+m a^T\), and

$$ E[b_j(X)b_j(X)^T] =G_jMG_j^T+G_jm g_j^T+g_jm^TG_j^T+g_jg_j^T. $$

The stochastic-integral terms have zero expectation: under the fourth-moment assumption, their integrands are square-integrable on finite time intervals, so the corresponding Itô integrals are true martingales. This justifies the expectation step rather than merely assuming it.

To obtain the covariance equation, subtract the derivative of \(mm^T\). Because \(m'=Am+a\),

$$ (mm^T)'=(Am+a)m^T+m(Am+a)^T. $$

Subtracting this from \(M'\) cancels the affine drift terms involving \(a\), leaving \(AP+PA^T\) from the drift. For the diffusion term, use \(M=P+mm^T\); then

$$ E[b_j(X)b_j(X)^T] =G_jPG_j^T+(G_jm+g_j)(G_jm+g_j)^T. $$

This proves both equations and shows why neither the diffusion term nor the distinction between \(M\) and \(P\) can be ignored. \(\square\)

Worked Checks

Worked Example: Additive Noise in a Mean-Reverting State

Let \(dX=-k(X-C)\,dt+\sigma\,dW\), with \(k>0\), and suppose \(m(0)=8\), \(P(0)=1\), \(C=10\), \(k=0.5\), and \(\sigma=2\). The mean equation is \(m'=-k(m-C)\), so

m(t)=C+(m(0)-C)e^{-kt} =10-2e^{-0.5t}.

For additive noise, the covariance equation is \(P'=-2kP+\sigma^2\). Solving this scalar linear equation gives

P(t)=e^{-2kt}P(0)+\frac{\sigma^2}{2k}(1-e^{-2kt}) =e^{-t}+4(1-e^{-t}).

At \(t=1\), \(m(1)=10-2e^{-0.5}\approx8.78694\), and \(P(1)=e^{-1}+4(1-e^{-1})\approx2.89636\). If the diffusion term is omitted, the erroneous equation is \(P'=-P\), which gives \(P(1)=e^{-1}\approx0.36788\). The discrepancy is not a small algebraic detail: it removes the variance continually injected by the noise. As a further check, the raw second moment is \(M(1)=P(1)+m(1)^2\approx80.10665\), not \(2.89636\).

Worked Example: Multiplicative Noise Is Not Additive Noise

Let \(dX=-kX\,dt+\gamma X\,dW\), with \(X(0)=5\), \(k=1\), and \(\gamma=0.8\). Here the diffusion coefficient is proportional to the current state. The mean equation gives \(m'=-m\), so \(m(t)=5e^{-t}\). The raw second-moment equation includes the square of the diffusion coefficient:

M'=(-2k+\gamma^2)M=(-2+0.64)M=-1.36M, \qquad M(t)=25e^{-1.36t}.

At \(t=1\), \(m(1)\approx1.83940\) and \(M(1)=25e^{-1.36}\approx6.41652\). Therefore

P(1)=M(1)-m(1)^2 \approx6.41652-(1.83940)^2 \approx3.0331.

A mistaken additive-noise calculation might freeze the diffusion at its initial value, \(\gamma X(0)=4\), and then use a constant diffusion coefficient of \(4\). That is a different SDE. Its variance equation would be \(P'=-2P+16\), giving \(P(1)=8(1-e^{-2})\approx6.91732\) from \(P(0)=0\). The value does not match the multiplicative-noise model. In that model the instantaneous diffusion changes with \(X(t)\), so the relevant term is \(\gamma^2E[X^2]\), not a fixed \(\sigma^2\).

Worked Example: A Two-Compartment Raw Moment Is Not a Covariance

Consider two states with independent noise:

dX_1=-X_1\,dt+dW_1,\qquad dX_2=-2X_2\,dt+2\,dW_2,\qquad (X_1(0),X_2(0))=(2,1).

The mean vector is \(m(t)=(2e^{-t},e^{-2t})^T\). Since the initial state is fixed and the noise drivers are independent, the covariance is diagonal. The scalar covariance equations are \(P_{11}'=-2P_{11}+1\) and \(P_{22}'=-4P_{22}+4\), both starting at zero; the cross-covariance remains zero. Hence

P(t)= \begin{pmatrix} \frac12(1-e^{-2t})&0\\ 0&1-e^{-4t} \end{pmatrix}.

At \(t=1\), this is approximately \(\begin{pmatrix}0.43233&0\\0&0.98168\end{pmatrix}\). But the raw second moment is \(M=P+mm^T\), so

M(1)= \begin{pmatrix} \frac12(1-e^{-2})+4e^{-2}&2e^{-3}\\ 2e^{-3}&1 \end{pmatrix} \approx \begin{pmatrix} 0.97367&0.09957\\ 0.09957&1.00000 \end{pmatrix}.

The nonzero raw cross-moment does not imply a nonzero covariance: here \(M_{12}=m_1m_2\), while \(P_{12}=M_{12}-m_1m_2=0\). Confusing these quantities can lead to a false conclusion about dependence between compartment states.

Diagnostics, Alternatives, and Limitations

A calculation is suspect if it gives zero noise contribution to a second moment when the diffusion is nonzero, or if it labels \(E[XX^T]\) as covariance without subtracting \(mm^T\). For multiplicative diffusion, check that the squared diffusion coefficient appears in the raw second-moment equation. In a scalar model \(dX=\alpha X\,dt+\gamma X\,dW\), for instance, \(M'=(2\alpha+\gamma^2)M\); a result containing \(2\alpha\) but not \(\gamma^2\) has lost quadratic variation.

Moment equations are exact under their stated assumptions and can be solved as deterministic ODEs. This is often more direct than estimating moments by Monte Carlo simulation, which introduces sampling error. Simulation remains useful as an independent diagnostic, but it does not repair a wrongly specified diffusion or an omitted term in the equations. Replacing multiplicative noise by additive noise can be a modelling approximation, not an algebraic simplification: it needs a stated approximation point and validation against the intended model.

These equations describe moments, not the full distribution of \(X(t)\). Even when the moment equations are closed, the distribution need not be summarized adequately by its mean and covariance. The derivation also relies on the affine coefficients and the stated integrability assumptions; nonlinear coefficients can require higher moments and need not produce a closed system.

Key takeaway: For linear SDEs, apply Itô's product rule before taking expectations, retain the quadratic-variation term, and recover covariance by subtracting the mean outer product from the raw second moment.

Check Your Understanding

Try each check, then compare with the solutions below.

For \(dX=-2X\,dt+3\,dW\), write the differential equation for \(P(t)\).
For \(dX=-X\,dt+0.5X\,dW\), with deterministic \(X(0)=x_0\), give the equations for \(m(t)\) and \(M(t)\).
If a calculation gives \(M=10\) and \(m=2\) for a scalar state, what is its variance?
In the two-compartment example, why can \(M_{12}\) be nonzero while \(P_{12}=0\)?

Solutions. (1) The additive-noise covariance equation is \(P'=-4P+9\). (2) The mean satisfies \(m'=-m\); the raw second moment satisfies \(M'=(-2+0.25)M=-1.75M\). (3) \(P=M-m^2=10-4=6\). (4) The raw cross-moment includes the product of the means: \(M_{12}=P_{12}+m_1m_2\). Independent noise and fixed initial states give \(P_{12}=0\), but both means are nonzero at finite time, so \(M_{12}=m_1m_2\) can be nonzero.