Kramers–Moyal expansion

In stochastic processes, Kramers–Moyal expansion refers to a Taylor series expansion of the master equation, named after Hans Kramers and José Enrique Moyal.

Start with the integro-differential master equation

${\frac {\partial p(x,t)}{\partial t}}=\int p(x,t|x_{0},t_{0})p(x_{0},t_{0})dx_{0}$

where $p(x,t|x_{0},t_{0})$ is the transition probability function, and $p(x,t)$ is the probability density at time $t$.

The Kramers–Moyal expansion transforms the above to an infinite order partial differential equation

$\partial _{t}p(x,t)=\sum _{n=1}^{\infty }(-\partial _{x})^{n}[D_{n}(x,t)p(x,t)]$

and also$$\partial _{t}p(x,t|x_{0},t_{0})=\sum _{n=1}^{\infty }(-\partial _{x})^{n}[D_{n}(x,t)p(x,t|x_{0},t_{0})]$$

where $D_{n}(x,t)$ are the Kramers–Moyal coefficients, defined by$$D_{n}(x,t)={\frac {1}{n!}}\lim _{\tau \to 0}{\frac {1}{\tau }}\mu _{n}(t|x,t-\tau )$$and $\mu _{n}$ are the central moment functions, defined by

$\mu _{n}(t'|x,t)=\int _{-\infty }^{\infty }(x'-x)^{n}p(x',t'\mid x,t)\ dx'.$

The Fokker–Planck equation is obtained by keeping only the first two terms of the series in which $\alpha _{1}$ is the drift and $\alpha _{2}$ is the diffusion coefficient.

Proof

In usual probability, where the probability density does not change, the moments of a probability density function determines the probability density itself by a Fourier transform (details may be found at the characteristic function page):$$p(x)={\frac {1}{2\pi }}\int e^{-ikx}{\tilde {p}}(k)dk=\sum _{n=0}^{\infty }{\frac {(-1)^{n}}{n!}}\delta ^{(n)}(x)\mu _{n}$$$${\tilde {p}}(k)=\int e^{ikx}p(x)dx=\sum _{n=0}^{\infty }{\frac {(ik)^{n}}{n!}}\mu _{n}$$Similarly, $$p(x,t|x_{0},t_{0})=\sum _{n=0}^{\infty }{\frac {(-1)^{n}}{n!}}\delta ^{(n)}(x-x_{0})\mu _{n}(t|x_{0},t_{0})$$ Now we need to integrate away the Dirac delta function. Fixing a small $\tau >0$, we have by the Chapman-Kolmogorov equation,$${\begin{aligned}p(x,t)&=\int p(x,t|x',t-\tau )p(x',t-\tau )dx'\\&=\sum _{n=0}^{\infty }{\frac {(-1)^{n}}{n!}}\int p(x',t-\tau )\delta ^{(n)}(x-x')\mu _{n}(t|x',t-\tau )dx'\\&=\sum _{n=0}^{\infty }{\frac {(-1)^{n}}{n!}}\partial _{x}^{n}(p(x,t-\tau )\mu _{n}(t|x,t-\tau ))\end{aligned}}$$The $n=0$ term is just $p(x,t-\tau )$, so taking derivative with respect to time,$$\partial _{t}p(x,t)=\lim _{\tau \to 0^{+}}{\frac {1}{\tau }}\sum _{n=1}^{\infty }{\frac {(-1)^{n}}{n!}}\partial _{x}^{n}(p(x,t-\tau )\mu _{n}(t|x,t-\tau ))=\sum _{n=1}^{\infty }(-\partial _{x})^{n}(p(x,t)D_{n}(x,t))$$

The same computation with $p(x,t|x_{0},t_{0})$ gives the other equation.

Forward and backward equations

The equation can be recast into a linear operator form, using the idea of infinitesimal generator.

Define the linear operator $${\mathcal {A}}f:=\sum _{n=1}^{\infty }(-\partial _{x})^{n}[D_{n}(x,t)f(x,t)]$$then the equation above states $${\begin{aligned}\partial _{t}p(x,t)&={\mathcal {A}}p(x,t)\\\partial _{t}p(x,t|x_{0},t_{0})&={\mathcal {A}}p(x,t|x_{0},t_{0})\end{aligned}}$$In this form, the equations are precisely in the form of a general Kolmogorov forward equation. The backward equation then states that$$\partial _{t}p(x_{1},t_{1}|x,t)=-{\mathcal {A}}^{\dagger }p(x_{1},t_{1}|x,t)$$where$${\mathcal {A}}^{\dagger }f:=\sum _{n=1}^{\infty }D_{n}(x,t)\partial _{x}^{n}[f(x,t)]$$

is the Hermitian adjoint of ${\mathcal {A}}$.

Computing the Kramers–Moyal coefficients

By definition,$$D_{n}(x,t)={\frac {1}{n!}}\lim _{\tau \to 0}{\frac {1}{\tau }}\mu _{n}(t|x,t-\tau )$$This definition works because $\mu _{n}(t|x,t)=0$, as those are the central moments of the Dirac delta function.

Since the even central moments are nonnegative, we have $D_{2n}\geq 0$ for all $n\geq 1$.

When the stochastic process is the Markov process $dX=bdt+\sigma dW_{t}$, we can directly solve for $p(x,t|x,t-\tau )$ as approximated by a normal distribution with mean $x+b(x)\tau$ and variance $\sigma ^{2}\tau$. This then allows us to compute the central moments, and so$$D_{1}=b,\quad D_{2}={\frac {1}{2}}\sigma ^{2},\quad D_{3}=D_{4}=\cdots =0$$This then gives us the 1-dimensional Fokker–Planck equation:$$\partial _{t}p=-\partial _{x}(bp)+{\frac {1}{2}}\partial _{x}^{2}(\sigma ^{2}p)$$

Pawula theorem

Let the operator $L_{m}$ be defined such $L_{m}f=\sum _{n=1}^{m}{\frac {(-1)^{n}}{n!}}{\frac {\partial ^{n}}{\partial x^{n}}}[\mu _{n}f]$. The probability density evolves by $\partial _{t}\rho \approx L_{m}\rho$. Different order of $m$ gives different level of approximation.

• $m=0$: the probability density does not evolve
• $m=1$: it evolves by deterministic drift only.
• $m=2$: it evolves by drift and Brownian motion (Fokker-Planck equation)
• $m=\infty$: the fully exact equation.

Pawula theorem states that for any other choice of $m$, there exists a probability density function $\rho$ that can become negative during its evolution $\partial _{t}\rho \approx L_{m}\rho$ (and thus fail to be a probability density function).

However, this doesn't mean that Kramers-Moyal expansions truncated at other choices of $m$ is useless. Though the solution must have negative values at least for sufficiently small times, the resulting approximation probability density may still be better than the $m=2$ approximation.

Proof

By Cauchy–Schwarz inequality, the central moment functions satisfy $\mu _{n+m}^{2}\leq \mu _{2n}\mu _{2m}$. So, taking the limit, we have $D_{n+m}^{2}\leq {\frac {(2n)!(2m)!}{(n+m)!^{2}}}D_{2n}D_{2m}$.

If some $D_{2+n}\neq 0$ for some $n\geq 1$, then $D_{2}D_{2+2n}>0$. In particular, $D_{2+n},D_{2+2n},D_{2+4n},...>0$. So the existence of any nonzero coefficient of order $\geq 3$ implies the existence of nonzero coefficients of arbitrarily large order.

Also, if $D_{n}\neq 0$, then $D_{2}D_{2n-2}>0,D_{4}D_{2n-4}>0,...$. So the existence of any nonzero coefficient of order $n$ implies all coefficients of order $2,4,...,2n-2$ are positive.

Thus, either the sequence $D_{1},D_{2},D_{3},...$ becomes zero at the third term, or all its even terms are positive.

Implementations

• Implementation as a python package