On the interpretation of Stratonovich calculus

Stochastic dynamical models are basic to the understanding of the role of random forcing in a wide range of scientific and engineering systems (e.g. [2–4]). The central theoretical approaches arise from the Einstein and Langevin studies of Brownian motion [5, 6], which provide mathematically different but physically equivalent and complimentary descriptions of the fate of a body (a pollen particle in water observed under a microscope) under the influence of random non-deterministic collisions (water molecules). Einstein determined the time evolution of the probability density of particles by solving the Fokker–Planck equation whereas Langevin wrote down an explicit deterministic momentum equation for a particle augmented by a Gaussian white noise forcing which perturbs the particle trajectory. The approach of Langevin now constitutes a canonical 'stochastic differential equation' (SDE) for a daunting scope of systems, but since its introduction the mathematical and physical interpretation of the noise term has been discussed and debated.

Whether viewing the problem from configuration space through the solution of the Fokker–Planck equation, solving a Langevin equation, or generating statistical realizations by performing Monte-Carlo simulations of the particles under consideration (e.g. [7]), a consistent interpretation and calculational scheme of the noise structure is needed. The approach developed by Itô [8] rests upon the Markovian and Martingale properties; the former captures the concept of a 'memoryless' process, wherein the conditional probability distribution of future states depends solely on the present state and the latter that, given all prior events, the expectation value of future stochastic events equals the present value. These properties have the advantage of simplifying many complicated time integrals but the disadvantage of requiring a new calculus which does not obey the traditional chain rule. In contrast, the approach of Stratonovich [9] does not invoke the Martingale property and preserves the chain rule and allows white noise to be treated as a regular derivative of a Brownian (or Weiner) process, ${{W}_{t}}$.

A model SDE (to which we shall return later) for a variable x(t) is

$$\begin{eqnarray}&&\frac{\text{d}x}{\text{d}t}=a\left( x,t \right)+b\left( x,t \right)\frac{\text{d}{{W}_{t}}}{\text{d}t},\end{eqnarray} \tag{ 1 }$$

wherein the first term is 'deterministic' and the second term is 'stochastic'. The case in which $b\left( x,t \right)$ is constant is referred to as additive noise and when it depends on $x=x\left( t \right)$ it is multiplicative noise. Now, despite being able to integrate equation (1) in a formal sense as

$$\begin{eqnarray}&&x\left( t \right)-x\left( {{t}_{0}} \right)=\mathop{\int }_{{{t}_{0}}}^{t}a\left( x,{{t}^{\prime }} \right)\text{d}{{t}^{\prime }}+\mathop{\int }_{{{t}_{0}}}^{t}b\left( x,{{t}^{\prime }} \right)\text{d}{{W}_{t}}\left( {{t}^{\prime }} \right),\end{eqnarray} \tag{ 2 }$$

the crux of the 'Itô–Stratonovich dilemma' resides in the integral over the Brownian process ${{W}_{t}}$. The issue is laid bare by recalling the definition of the definite integral of a real valued function f(t) in terms of the Riemann sum viz.,

$$\begin{eqnarray}&&\mathop{\int }_{{{t}_{0}}}^{t}f\left( {{t}^{\prime }} \right)\text{d}{{t}^{\prime }}\equiv \underset{n\to \infty }{\mathop{\lim }}\,\underset{k=0}{\overset{n-1}{\mathop{\mathop{\sum }}}}\,f\left( {{\tilde{t}\,}_{k}} \right)\left( {{t}_{k+1}}-{{t}_{k}} \right),\end{eqnarray} \tag{ 3 }$$

where importantly ${{\tilde{t}\,}_{k}}\in \left[ {{t}_{k}},{{t}_{k+1}} \right]$. The issues are (i) the Brownian process ${{W}_{t}}$ is nowhere differentiable and because $\frac{d{{W}_{t}}}{dt}$ is δ-autocorrelated then in any interval on the real line the white noise it represents fluctuates an infinite number of times with an infinite variance, and (ii) the defining limit of the integral in equation (3) depends on the place in $\left[ {{t}_{k}},{{t}_{k+1}} \right]$ where ${{\tilde{t}\,}_{k}}$ is chosen; the choice that provides the seed of the dilemma and the origin of the two different calculi.

In the Itô approach the choice is ${{\tilde{t}\,}_{k}}={{t}_{k}}$, which maintains the Martingale property due to the fact that ${{t}_{k}}$ is the present value in the integrand, thereby forcing the expectation value in $\left[ {{t}_{k}},{{t}_{k+1}} \right]$ to take the present value. In contrast, the Stratonovich approach defines ${{\tilde{t}\,}_{k}}=\left( {{t}_{k}}+{{t}_{k+1}} \right)/2$, which abandons the Martingale property and hence maintains the normal rules of calculus. Stratonovich referred to his choice as a 'symmetrization' between past and future. Despite its wide usage, we have not found a physical interpretation of the basis of the past/future symmetrization of Stratonovich calculus.

The outline of this paper is as follows. In the next section we summarize our approach and discuss its connection to the dilemma generally and to other studies. We then outline the conventional physical science viewpoint and the perspective of stochastic calculus, before coming to our main point. Finally, we discuss the relationship to other approaches before concluding.

Wong and Zakai [1] argued that in any real world system perfect white noise does not exist and that Brownian motion x(t) approximates a description ${{x}_{n}}\left( t \right)$ that is continuous with at least a piece-wise continuous derivative. By showing that ${{x}_{n}}\left( t \right)\to x\left( t \right)$ as $n\to \infty$ they recovered Stratonovich calculus. Accordingly, the choice of stochastic calculus resides in the characteristics of the noise and continuity arguments. In finance, it is argued that short time scale processes are truly discontinuous and thus Itô calculus is preferred (e.g. [10]), thereby maintaining the Martingale property. However, in physics the continuous motion of Brownian particles influenced by high frequency white noise has long been considered within the framework of normal calculus [11]. Conceptually, there is no clear distinction between the statistics of water molecules colliding pollen grains and trading options or stocks. Hence, it remains an open exercise to answer the questions of when, how and whether to use continuity considerations as a core criterion to choose either of the calculi being discussed here.

As noted in the introduction, the complimentary approaches of Einstein and Langevin, with a reliance on δ-autocorrelated noise, provide consistent testable predictions of Brownian motion. The combination of this and the additional consistency with Stratonovich calculus led to the suggestion that physical scientists should avoid Itô calculus [12]. Although Van Kampen [13] cautiously suggested that the Langevin equation is intrinsically insufficient for representing systems with internal noise, it is still the case that Stratonovich calculus is appropriate for both internal and external noise and there are an enormous number of problems in which one cannot make this distinction.

Here we focus on an ostensibly physical argument to discuss the origin of Stratonovich calculus. We generalize the theorem Wong and Zakai [1] using the integral Taylor expansion and appealing solely to the ${{L}^{2}}$ integrability of a function, thereby avoiding discussions of the regularity of a stochastic process. Because the principal difference between the Itô–Langevin and the Stratonovich–Langevin formulations lies in the drift term, we focus on this in our approach and examine the intrinsic nature of short time scale processes approximated by white noise to decide which calculus is more appropriate. As opposed to Turelli [14], who took as a continuous deterministic model an approximation of a discretized model (in population dynamics), we begin with a well defined deterministic model.

The term 'white noise' refers to a random signal $\Gamma \left( t \right)$ in which all frequencies contribute equally to the power spectral density. For a general stationary random process the autocorrelation $\mathcal{R}\left( {{t}^{\prime }}-t \right)=\langle\Gamma \left( {{t}^{\prime }} \right)\Gamma \left( t \right)\rangle$, where $\langle\cdot \rangle$ is the ensemble average, is by definition independent of time and symmetric/even $\mathcal{R}\left( {{t}^{\prime }}-t \right)=\mathcal{R}\left( t-{{t}^{\prime }} \right)$. If $t-{{t}^{\prime }}\equiv\tau$, then let ${{\tau }_{c}}$ be the de-correlation time such that $\mathcal{R}\left( {{\tau }_{c}} \right)={{\text{e}}^{-1}}\mathcal{R}\left( 0 \right)$. Typically, as ${{\tau }_{c}}\to 0$ and $\mathop{\int }_{-\infty }^{\infty }\mathcal{R}\left( {{\tau }_{c}} \right)\text{d}{{\tau }_{c}}$ remains finite, we say that $\mathcal{R}\left( {{\tau }_{c}} \right)=\delta \left( 0 \right)$ and the white noise has a δ-autocorrelation structure.

In science and engineering the traditional form of the Langevin equation, slightly different than that in equation (1), is written as

$$\begin{eqnarray}&&\frac{\text{d}x}{\text{d}t}=a\left( x,t \right)+b\left( x,t \right)\Gamma \left( t \right),\end{eqnarray} \tag{ 4 }$$

where the noise forcing is δ-autocorrelated as $\langle\Gamma \left( {{t}^{\prime }} \right)\Gamma \left( t \right\rangle=\delta \left( {{t}^{\prime }}-t \right)$. Based on this, we calculate $\langle\frac{\text{d}x}{\text{d}t}\rangle$ starting with the first integral

$$\begin{eqnarray}&&x\left( t+\tau \right)-x\left( t \right)=\mathop{\int }_{t}^{t+\tau }\left( a\left[ x\left( {{t}^{\prime }} \right),{{t}^{\prime }} \right]+b\left[ x\left( {{t}^{\prime }} \right),{{t}^{\prime }} \right]\Gamma \left( {{t}^{\prime }} \right) \right)\text{d}{{t}^{\prime }}.\end{eqnarray} \tag{ 5 }$$

We assume that τ is small relative to the time over which the macroscopic system evolves, but large relative to the time scale of fluctuations in the stochastic process $\Gamma \left( t \right)$. When $a\left( x,t \right)$ and $b\left( x,t \right)$ are analytic, we can Taylor expand them about a time t as

$$\begin{eqnarray}\begin{matrix} a\left( x\left( {{t}^{\prime }} \right),{{t}^{\prime }} \right) & = & a\left( x\left( t \right),t \right)+\frac{\partial a}{\partial x}\left( x\left( {{t}^{\prime }} \right)-x\left( t \right) \right)+\frac{\partial a}{\partial t}\left( {{t}^{\prime }}-t \right)+\ldots \\ b\left( x\left( {{t}^{\prime }} \right),{{t}^{\prime }} \right) & = & b\left( x\left( t \right),t \right)+\frac{\partial b}{\partial x}\left( x\left( {{t}^{\prime }} \right)-x\left( t \right) \right)+\frac{\partial b}{\partial t}\left( {{t}^{\prime }}-t \right)+\ldots . \\ &&\end{matrix}\end{eqnarray} \tag{ 6 }$$

following Risken [15]. Inserting (6) into (5) leads to

$$\begin{eqnarray}\begin{matrix} x\left( t+\tau \right)-x\left( t \right) & = & a\left( x\left( t \right),t \right)\tau +\frac{\partial a}{\partial x}\mathop{\int }_{t}^{t+\tau }\left( x\left( {{t}^{\prime }} \right)-x\left( t \right) \right)\text{d}{{t}^{\prime }}+\frac{\partial a}{\partial t}\frac{1}{2}{{\tau }^{2}} \\ {} & {} & +b\left( x\left( t \right),t \right)\mathop{\int }_{t}^{t+\tau }\Gamma \left( {{t}^{\prime }} \right)\text{d}{{t}^{\prime }}+\frac{\partial b}{\partial x}\mathop{\int }_{t}^{t+\tau }\left( x\left( {{t}^{\prime }} \right)-x\left( t \right) \right)\Gamma \left( {{t}^{\prime }} \right)\text{d}{{t}^{\prime }} \\ {} & {} & +\frac{\partial b}{\partial t}\mathop{\int }_{t}^{t+\tau }\left( {{t}^{\prime }}-t \right)\Gamma \left( {{t}^{\prime }} \right)\text{d}{{t}^{\prime }}+h.o.t. \\ &&\end{matrix}\end{eqnarray} \tag{ 7 }$$

For $x\left( {{t}^{\prime }} \right)-x\left( t \right)$, we can treat equation (5) iteratively and then take the ensemble average, which is

$$\begin{eqnarray}&&\langle x\left( t+\tau \right)-x\left( t \right)\rangle =a\left( x,t \right)\tau +b\frac{\partial b}{\partial x}\mathop{\int }_{t}^{t+\tau }\;\text{d}{{t}^{\prime }}\mathop{\int }_{t}^{{{t}^{\prime }}}\langle \Gamma \left( {{t}^{\prime\prime }} \right)\Gamma \left( {{t}^{\prime }} \right)\rangle \text{d}{{t}^{\prime\prime }}+O\left( {{\tau }^{2}} \right).\end{eqnarray} \tag{ 8 }$$

Therefore,

$$\begin{eqnarray}&&\underset{\tau \to 0}{\mathop{\lim }}\,\frac{1}{\tau }\langle x\left( t+\tau \right)-x\left( t \right)\rangle =a\left( x,t \right)+b\frac{\partial b}{\partial x}\frac{1}{\tau }\mathop{\int }_{t}^{t+\tau }\;\text{d}{{t}^{\prime }}\mathop{\int }_{t}^{{{t}^{\prime }}}\langle \Gamma \left( {{t}^{\prime\prime }} \right)\Gamma \left( {{t}^{\prime }} \right)\rangle \text{d}{{t}^{\prime\prime }}.\end{eqnarray} \tag{ 9 }$$

By the definition of white noise discussed above, $\langle \Gamma \left( {{t}^{\prime }} \right)\Gamma \left( {{t}^{\prime\prime }} \right)\rangle =\delta \left( {{t}^{\prime }}-{{t}^{\prime\prime }} \right)$, we have

$$\begin{eqnarray}&&\mathop{\int }_{t}^{t+\tau }\;\text{d}{{t}^{\prime }}\mathop{\int }_{t}^{{{t}^{\prime }}}\langle \Gamma \left( {{t}^{\prime\prime }} \right)\Gamma \left( {{t}^{\prime }} \right)\rangle \text{d}{{t}^{\prime\prime }}=\mathop{\int }_{t}^{t+\tau }\;\text{d}{{t}^{\prime }}\mathop{\int }_{t}^{{{t}^{\prime }}}\delta \left( {{t}^{\prime\prime }}-{{t}^{\prime }} \right)\text{d}{{t}^{\prime\prime }}.\end{eqnarray} \tag{ 10 }$$

In the integral $\mathop{\int }_{t}^{{{t}^{\prime }}}\delta \left( {{t}^{\prime\prime }}-{{t}^{\prime }} \right)\text{d}{{t}^{\prime\prime }}$, half of the contribution of the delta function is considered such that the value of the integral becomes $\frac{1}{2}$. Therefore, $\mathop{\int }_{t}^{t+\tau }\;\text{d}{{t}^{\prime }}\mathop{\int }_{t}^{{{t}^{\prime }}}\langle \Gamma \left( {{t}^{\prime\prime }} \right)\Gamma \left( {{t}^{\prime }} \right)\rangle \text{d}{{t}^{\prime\prime }}$ is equal to $\frac{1}{2}\tau$. What is the physical meaning of the half contribution of the delta function? The entire interval $\left[ t,{{t}^{\prime }} \right]$ does not contribute to the integral $\mathop{\int }_{t}^{{{t}^{\prime }}}\delta \left( {{t}^{\prime\prime }}-{{t}^{\prime }} \right)\text{d}{{t}^{\prime\prime }}$. We shrink the interval near ${{t}^{\prime }}$ as $\left[ {{t}^{\prime }}-\epsilon ,{{t}^{\prime }} \right]$, where can be interpreted as the size of a subinterval in the Riemann summation discussed above in distinguishing the stochastic calculi. Because the delta function is even $\delta \left( {{t}^{\prime\prime }}-{{t}^{\prime }} \right)=\delta \left( {{t}^{\prime }}-{{t}^{\prime\prime }} \right)$, the half contribution can be written as

$$\begin{eqnarray}&&\mathop{\int }_{t}^{{{t}^{\prime }}}\delta \left( {{t}^{\prime\prime }}-{{t}^{\prime }} \right)\text{d}{{t}^{\prime\prime }}=\frac{1}{2}\mathop{\int }_{{{t}^{\prime }}-\epsilon }^{{{t}^{\prime }}+\epsilon }\delta \left( {{t}^{\prime\prime }}-{{t}^{\prime }} \right)\text{d}{{t}^{\prime\prime }}.\end{eqnarray} \tag{ 11 }$$

Despite this being just a simple mathematical modification, we can interpret it as a representation that the weighting from past is equal to that from future in the neighbourhood of ${{t}^{\prime }}$. The origin of this symmetry comes from the time symmetry of the autocorrelation function $\mathcal{R}\left( \tau \right)$. We will connect this half contribution of the autocorrelation to the Stratonovich integral below, but we note now the above development yields

$$\begin{eqnarray}&&\left\langle \frac{\text{d}x}{\text{d}t} \right\rangle =a\left( x,t \right)+\frac{1}{2}b\left( x,t \right)\frac{\partial }{\partial x}b\left( x,t \right),\end{eqnarray} \tag{ 12 }$$

which is exactly same as that obtained from Stratonovich calculus.

Mathematicians have long questioned the validity of the conventional definition of white noise because Brownian motion ${{W}_{t}}$ is nowhere differentiable. Namely, defining $\Gamma \left( t \right)$ as $\frac{\text{d}{{W}_{t}}}{\text{d}t}$ is poorly grounded in normal calculus and we must revisit the interpretation of equation (4) in light of the mathematical definition of Brownian motion ${{W}_{t}}$;

• (i)  
  ${{W}_{0}}$ = 0
• (ii)  
  ${{W}_{t}}$ is almost surely continuous (or sample-continuous) with respect to t
• (iii)  
  ${{W}_{t}}$ has independent increments and ${{W}_{t}}-{{W}_{s}}$, with $s<t$, is governed by a normal distribution, $N\left( 0,t-s \right)$, with zero mean and standard deviation $\sqrt{t-s}$.

Therefore, equation (4) can be rewritten as

$$\begin{eqnarray}&&\text{d}x=a\left( x,t \right)\text{d}t+b\left( x,t \right)\text{d}{{W}_{t}}.\end{eqnarray} \tag{ 13 }$$

On the one hand, using this integral form avoids the issue of the differentiability of ${{W}_{t}}$, but this comes at the cost of having to introduce a different sort of calculus. However, before coming to the issue of the implementation of a particular stochastic calculus, we perform the same expansion as discussed previously, in (7) and (8), which leads to

$$\begin{eqnarray}&&\langle x\left( t+\tau \right)-x\left( t \right)\rangle =a\left( x,t \right)\tau +b\left( x,t \right)\frac{\partial }{\partial x}b\left( x,t \right)\left\langle \mathop{\int }_{t}^{t+\tau }\mathop{\int }_{t}^{{{t}^{\prime }}}\text{d}W_{t}^{\prime\prime }\text{d}W_{t}^{\prime } \right\rangle +O\left( {{\tau }^{2}} \right).\end{eqnarray} \tag{ 14 }$$

The first term $a\left( x,t \right)$ can be interpreted as the deterministic response to macro-scale forcing and the second term embodies the interaction between a macroscopic quantity $b\left( x,t \right)$ and the microscopic accumulation of random noise forcing. On a time τ that is short from the macroscopic perspective, the interaction between two scales manifests itself multiplicatively. The contribution of the micro-scale to the change of the macroscopic quantity x appears as a double integral of the Brownian motion. It is at this point that we must decide upon a specific calculus, and this choice determines the value of the double integral in equation (14).

The double integral can be written as $\mathop{\int }_{t}^{t+\tau }{{W}^{\prime }}\text{d}{{W}^{\prime }}$ and approximated by a Riemann sum as

$$\begin{eqnarray}&&\mathop{\int }_{t}^{t+\tau }{{W}^{\prime }}\text{d}{{W}^{\prime }}=\underset{k=0}{\overset{N-1}{\mathop{\mathop{\sum }}}}\,{{W}_{t_{k}^{\prime }}}\Delta {{W}_{{{t}_{k}}}},\end{eqnarray} \tag{ 15 }$$

where $\Delta {{W}_{{{t}_{k}}}}\equiv {{W}_{{{t}_{k+1}}}}-{{W}_{t_{k}^{\prime }}}$. As should be intuitive from the previous discussion, the confusion lies in the choice of $t_{k}^{\prime }$ for ${{W}_{t_{k}^{\prime }}}$. As can be seen below, this choice determines the average of the integral viz.,

$$\begin{eqnarray}\begin{matrix} \left\langle \underset{k=0}{\overset{k=N-1}{\mathop{\mathop{\sum }}}}\,{{W}_{t_{k}^{\prime }}}\Delta {{W}_{{{t}_{k}}}} \right\rangle & = & \underset{k=0}{\overset{k=N-1}{\mathop{\mathop{\sum }}}}\,\left\langle \left( {{W}_{t_{k}^{\prime }}}-{{W}_{{{t}_{k}}}}+{{W}_{{{t}_{k}}}} \right)\left( {{W}_{{{t}_{k+1}}}}-{{W}_{{{t}_{k}}}} \right) \right\rangle \\ {} & = & \underset{k=0}{\overset{k=N-1}{\mathop{\mathop{\sum }}}}\,\left\langle \left( {{W}_{t_{k}^{\prime }}}-{{W}_{{{t}_{k}}}} \right)\left( {{W}_{{{t}_{k+1}}}}-{{W}_{{{t}_{k}}}} \right) \right\rangle \\ {} & = & \underset{k=0}{\overset{k=N-1}{\mathop{\mathop{\sum }}}}\,\left( t_{k}^{\prime }-{{t}_{k}} \right). \\ &&\end{matrix}\end{eqnarray} \tag{ 16 }$$

Itô [8] argued that by choosing $t_{k}^{\prime }$ = ${{t}_{k}}$ the most important characteristic preserved is the Martingale property. This implies that the expectation value of a future event is equivalent to that of the present value in the subinterval $\left[ {{t}_{k}},{{t}_{k+1}} \right]$. Hence, the expectation value of the double integral is zero. The disadvantage of Itô calculus is that the chain rule of normal calculus is lost; one needs additional terms to preserve the Martingale property during integration. On the other hand, Stratonovich [9] chose $t_{k}^{\prime }$ = $\frac{{{t}_{k}}+{{t}_{k+1}}}{2}$, and the integral has a value of $\frac{1}{2}\tau$. This approach preserves the chain rule, implying that Stratonovich calculus must be used when we approximate noise processes that are continuous and have piece-wise continuous derivatives (e.g. [1]). When we approximate the rapid fluctuations of a stochastic process as white noise, then we must preserve the normal rules of calculus if we are to capture the dynamics of Brownian motion.

In a classical physical world where the continuity is more or less guaranteed, Stratonovich calculus is evidently the most appropriate. However, when we consider the colliding of particles or the buying or selling of stocks or options, these short time-scale phenomena appear to be discontinuous. Such circumstances lead us to question the validity of the continuity of a stochastic process. That said, it may be impossible to find events that have zero autocorrelation. Thus it is a natural question to examine the effect of a non-zero autocorrelation in a multiplicative noise process by investigating limiting behaviors, which may provide some insight into the nature of the appropriate stochastic calculus. Hence we consider the simplest case of multiplicitive behavior; $\text{d}x/\text{d}t=xF\left( t \right)$, where F(t) is a stochastic process with a finite decorrelation time. For this example, using the same procedure used to arrive at equation (8), we find that

$$\begin{eqnarray}&&\langle x\left( t+\frac{\tau }{n} \right)-x\left( t \right)\rangle =x\left( t \right)\mathop{\int }_{t}^{t+\frac{\tau }{n}}\mathop{\int }_{t}^{{{t}^{\prime }}}F\left( {{t}^{\prime\prime }} \right)\langle F\left( {{t}^{\prime }} \right)\text{d}{{t}^{\prime\prime }}\rangle \text{d}{{t}^{\prime }}.\end{eqnarray} \tag{ 17 }$$

An important feature of our argument involves the time $\tau /n$, with integer n, which is a subdivision of the total domain $\left[ t,t+\tau \right]$. For example, consider the commonly used autocorrelation $\langle F\left( {{t}^{\prime\prime }} \right)F\left( {{t}^{\prime }} \right)\rangle =\frac{{{\sigma }^{2}}}{2{{\tau }_{c}}}{{\text{e}}^{-|{{t}^{\prime\prime }}-{{t}^{\prime }}|/{{\tau }_{c}}}}$, where ${{\tau }_{c}}$ is the decorrelation time of the amplitude of the stochastic noise F(t) satisfying ${{\tau }_{c}}\tau$. Now, let ${{\tau }_{\text{m}}}\equiv \tau /n$ be a minimum time over which we can observe the decay of the stochastic noise.

Whereas it is always the case that $\tau \gg {{\tau }_{c}}$, the value of ${{\tau }_{\text{m}}}$ depends on the fineness n of the subdivision, with ${{\tau }_{\text{m}}}$ increasing as n decreases; coarsening the temporal resolution. So long as ${{\tau }_{\text{m}}}\gg {{\tau }_{c}}$ we can still justify the use of the stochastic differential equation but once ${{\tau }_{\text{m}}}\leqslant {{\tau }_{c}}$ this is no longer the case. Therefore, increasing ${{\tau }_{\text{m}}}$ to values above ${{\tau }_{c}}$ but still small relative to τ one reaches a value of ${{\tau }_{\text{m}}}\equiv {{\tau }_{\star }}$ that is sufficiently short that the decay of the noise can be observed but sufficiently large so that the stochastic differential equation is valid. This is achieved as follows. Consider n such that $\tau\gg {{\tau }_{\text{m}}}={{\tau }_{c}}$, which we rewrite as $1\gg\frac{1}{n}=\frac{{{\tau }_{c}}}{\tau }$. Now, as n decreases from $\frac{\tau }{{{\tau }_{c}}}$ it will reach a value $n={{n}_{\star }}$ such that $1\gg\frac{1}{{{n}_{\star }}}\gg\frac{{{\tau }_{c}}}{\tau }$ and our analysis is valid for ${{n}_{\star }}\geqslant n\gg 1$ as described presently.

With the above form of the autocorrelation, equation (17) becomes

$$\begin{eqnarray}&&\langle x\left( t+\frac{\tau }{n} \right)-x\left( t \right)\rangle =x\left( t \right)\times \frac{1}{2}{{\sigma }^{2}}\frac{\tau }{n}\left[ 1+\frac{{{\tau }_{c}}}{\tau /n}\left( {{e}^{-\frac{\tau }{n{{\tau }_{c}}}}}-1 \right) \right],\end{eqnarray} \tag{ 18 }$$

which for ${{n}_{\star }}\geqslant n\gg 1$ can be approximated as $\langle x\left( t+\frac{\tau }{n} \right)\rangle =x\left( t \right)\left( 1+\frac{1}{2}{{\sigma }^{2}}\frac{\tau }{n} \right)$ and then rewritten as

$$\begin{eqnarray}&&\langle x\left( t+\tau \right)\rangle =x\left( t \right){{\left( 1+\frac{1}{2}{{\sigma }^{2}}\frac{\tau }{n} \right)}^{n}}.\end{eqnarray} \tag{ 19 }$$

The limit ${{n}_{\star }}\geqslant n\gg 1$ insures that the right hand side converges to $x\left( t \right){{e}^{\frac{1}{2}{{\sigma }^{2}}\tau }}$ and hence

$$\begin{eqnarray}&&\langle x\left( t+\tau \right)\rangle =x\left( t \right){{e}^{\frac{1}{2}{{\sigma }^{2}}\tau }}\approx x\left( t \right)\left( 1+\frac{1}{2}{{\sigma }^{2}}\tau \right),\end{eqnarray} \tag{ 20 }$$

consistent with the result from Stratonovich calculus, in which $\frac{1}{2}{{\sigma }^{2}}\tau$ is analogous to $\langle {{W}_{t}}\Delta {{W}_{t}}\rangle$ in the Riemann sum with ${{W}_{t}}={{W}_{\left( {{t}_{k}}+{{t}_{k+1}} \right)/2}}$ and $\Delta {{W}_{t}}={{W}_{{{t}_{k+1}}}}-{{W}_{{{t}_{k}}}}$. This is the essence of Stratonovich calculus; choosing the mid-point in a subinterval is equivalent to including only half of the effect of a finite autocorrelation. Clearly, in this example, the magnitude of ${{\tau }_{c}}$ does not play a role unless ${{\tau }_{\text{m}}}\searrow {{\tau }_{c}}$ and the neccessary separation of time scales starts to be violated. Therefore, in this case the commonly-used definition of white noise $\langle F\left( {{t}^{\prime\prime }} \right)F\left( {{t}^{\prime }} \right)\rangle \propto \delta \left( {{t}^{\prime\prime }}-{{t}^{\prime }} \right)$ becomes an increasingly better approximation, which we can understand as follows.

Recall that

$$\begin{eqnarray}&&\mathop{\int }_{t}^{{{t}^{\prime }}}\frac{{{\sigma }^{2}}}{2{{\tau }_{c}}}{{e}^{-\frac{|{{t}^{\prime }}-{{t}^{\prime }}|}{{{\tau }_{c}}}}}\text{d}{{t}^{\prime }}\approx \mathop{\int }_{t}^{{{t}^{\prime }}}{{\sigma }^{2}}\delta \left( {{t}^{\prime\prime }}-{{t}^{\prime }} \right)\text{d}{{t}^{\prime\prime }}\sim \mathcal{O}\left( {{\sigma }^{2}} \right).\end{eqnarray} \tag{ 21 }$$

We can proceed along the same lines as above but without the need to invoke a specific form of the autocorrelation. Rather, we only need to assume that F(t) is stationary,

$$\begin{eqnarray}\begin{matrix} \mathop{\int }_{t}^{t+\frac{\tau }{n}}\mathop{\int }_{t}^{{{t}^{\prime }}}\langle F\left( {{t}^{\prime\prime }} \right)F\left( {{t}^{\prime }} \right)\rangle \text{d}{{t}^{\prime\prime }}\text{d}{{t}^{\prime }} & = & \mathop{\int }_{-\tau /n}^{0}\text{d}s\mathop{\int }_{t-s}^{t+\tau /n}K\left( s \right)\text{d}{{t}^{\prime }}=\mathop{\int }_{-\tau /n}^{0}K\left( s \right)\left( \tau +s \right)\text{d}s \\ {} & {} & \approx \frac{1}{2}\frac{\tau }{n}\mathop{\int }_{-\infty }^{\infty }K\left( s \right)\text{d}s, \\ &&\end{matrix}\end{eqnarray} \tag{ 22 }$$

where $K\left( s \right)=\langle F\left( t \right)F\left( t+s \right)\rangle$ and we exploited $\tau /n\gg {{\tau }_{c}}$, thereby neglecting s compared to $\tau /n$ (see e.g. chapter 15 of [16]). Here, $\frac{1}{2}\frac{\tau }{n}\mathop{\int }_{-\infty }^{\infty }K\left( s \right)\text{d}s$ is an approximation of $\langle {{W}_{t}}\Delta {{W}_{t}}\rangle$ in the Riemann sum, where the appropriate form of stochastic calculus is determined by the magnitude of $\mathop{\int }_{-\infty }^{\infty }K\left( s \right)\text{d}s$. Namely, if

$$\begin{eqnarray}&&\mathop{\int }_{-\infty }^{\infty }K\left( s \right)\text{d}s\sim \mathcal{O}\left( {{\sigma }^{2}} \right),\end{eqnarray} \tag{ 23 }$$

where $\langle {{F}^{2}}\left( t \right)\rangle ={{\sigma }^{2}}$, the determination of $\langle {{W}_{t}}\Delta {{W}_{t}}\rangle$ follows Stratonovich calculus and hence choosing a δ-function for K(s) is appropriate. However, if

$$\begin{eqnarray}&&\mathop{\int }_{-\infty }^{\infty }K\left( s \right)\text{d}s{{\sigma }^{2}},\end{eqnarray} \tag{ 24 }$$

then $\frac{1}{2}\frac{\tau }{n}\mathop{\int }_{-\infty }^{\infty }K\left( s \right)\text{d}s$ is negligible and the determination of $\langle {{W}_{t}}\Delta {{W}_{t}}\rangle$ follows Itô-calculus. Here, K(s) can be considered approximately as a zero-measure function in ${{L}^{2}}$ space such that a δ-function cannot be used for the approximation of the autocorrelation function K(s). This argument is similar to Moritaʼs [17], who also considered three time scales and his choice of stochastic calculus also relies on the condition for the autocorrelation function, but as manifests itself in the form of a generalized Fokker–Planck equation.

The choice of Stratonovich calculus for colored noise has been considered within the framework of a master equation by Horsthemke and Lefever [18]. Colored noise converges to Gaussian white noise through the application of the central limit theorem, which they argued through a rescaling of time ($\Delta t$) and space ($\Delta x$) assuming a diffusive process, viz., $\Delta t\to \Delta t/{{\epsilon }^{2}},\Delta x\to \Delta x/\epsilon$. In the limit $\epsilon \to 0$ the master equation becomes the Fokker–Planck equation with Gaussian white noise interpreted using Stratonovich calculus. Horsthemke and Lefever note that the interpretation of Stratonovich calculus in the approach of Wong and Zakai [1] requires the strong regularity of colored noise (a ${{\text{C}}_{1}}$ function) for convergence to Gaussian white noise. However, Wong and Zakai [1] demonstrate the convergence directly from an SDE, whereas Horsthemke and Lefever [18] avoid the issue of regularity by beginning with an Uhlenbeck–Ornstein process, which is not C$_{1}$, and demonstrate convergence using a master equation, but do not connect the finite autocorrelation of colored noise with a physical interpretation that underlies the construction of Stratonovich calculus.

It is the mid-point selection rule in the Riemann sum that is associated with the finite noise autocorrelation. However, traditionally this selection rule has been interpreted as the 'magical' choice required in order to recover the normal rules of calculus, but Stratonovich calculus itself does not provide any new dynamical insight that plays a role analogous to the Martingale property of Itô calculus. Hence, we seek a basic understanding of how the mid-point selection rule emerges as a consequence of applying the central limit theorem to colored noise. Our approach is to consider shot noise, which has an autocorrelation represented by a δ-function.

What is the origin of the δ-function? Brownian motion can be generated by a collection of independent random processes. In the Riemann-summation approximation of the integral $\mathop{\int }_{t}^{t+\tau }{{W}^{\prime }}\text{d}{{W}^{\prime }}$, we first considered the subdivision of the time domain $\left[ t,t+\tau \right]$ using $\tau /n$ with integer n. As discussed in the argument leading to equation (20), so long as it is large the freedom to choose n determined the time scale $\tau /n$ over which the decay of the noise forcing is observed; when $\tau /n\gg {{\tau }_{c}}$ then the noise signal appears as a discrete packet, appropriately simulated as shot noise (e.g. [19]).

Consider two 'square' signals with amplitude h and time duration w, occurring with probability $\alpha w$. When both are contained in the same packet and have the same sign ($\pm h$) they are positively correlated. Otherwise they are independent and uncorrelated. In the probability domain $\Omega$ the autocorrelation $\mathcal{R}\left( \tau \right)$ is

$$\begin{eqnarray}&&\mathcal{R}\left( \tau \right)={{\mathop{\int }}_{\Omega }}{{n}_{1}}{{n}_{2}}P\left( {{n}_{1}},{{n}_{2}} \right)\text{d}{{n}_{1}}\text{d}{{n}_{2}},\end{eqnarray} \tag{ 25 }$$

where ${{n}_{1}}$ (${{n}_{2}}$) is a random variable occurring at the time t ($t+\tau$) and $P\left( {{n}_{1}},{{n}_{2}} \right)$ is the joint probability density function of ${{n}_{1}}$ and ${{n}_{2}}$. We need only consider the case when τ is smaller than w, which is

$$\begin{eqnarray}\begin{matrix} \mathcal{R}\left( \tau \right) & = & {{\mathop{\int }}_{\Omega }}n_{1}^{2}P\left( {{n}_{1}} \right)P\left( {{n}_{2}}|{{n}_{1}} \right)\text{d}{{n}_{1}}\text{d}{{n}_{2}} \\ {} & = & \left( 1-\frac{\tau }{w} \right)\alpha w\left( \frac{1}{2}{{\left( \frac{h}{w} \right)}^{2}}+\frac{1}{2}{{\left( -\frac{h}{w} \right)}^{2}} \right) \\ {} & = & \alpha {{h}^{2}}\frac{1}{w}\left( 1-\frac{\tau }{w} \right), \\ &&\end{matrix}\end{eqnarray} \tag{ 26 }$$

from which we can readily calculate

$$\begin{eqnarray}&&\mathop{\int }_{-\infty }^{\infty }R\left( \tau \right)\text{d}\tau \approx \alpha {{h}^{2}}\mathop{\int }_{-w}^{w}\frac{1}{w}\left( 1-\frac{\tau }{w} \right)\text{d}\tau =2\alpha {{h}^{2}}.\end{eqnarray} \tag{ 27 }$$

Thus, the constant $\alpha {{h}^{2}}$ corresponds to ${{\sigma }^{2}}$, $h/w$ is analogous to $\text{d}W/\text{d}t$, and hence $\mathcal{R}\left( \tau \right)=2{{\sigma }^{2}}\delta \left( \tau \right)$. Therefore, despite our argument emerging from shot noise, it leads naturally to the definition of white noise in which there exists a finite time correlation that is too minute to be realized when viewed from a coarse time-scale, a coarseness which depends on the observer.

Clearly, a similar procedure applies directly to the Riemann-sum so long as the same constraints on the subdivision of the time domain are in place. When $\tau /n\simeq w$, two signals of the same sign are positively correlated in the interval. Hence, the random variables ${{n}_{1}}{{n}_{2}}$ only have the value ${{h}^{2}}$ with probability ${{\mathop{\int }}_{\Omega }}P\left( {{n}_{1}} \right)P\left( {{n}_{2}}|{{n}_{1}} \right)\text{d}{{n}_{1}}\text{d}{{n}_{2}}$, where $P\left( {{n}_{1}} \right)$ is $\alpha \tau /n$ and $P\left( {{n}_{2}}|{{n}_{1}} \right)$ is $1/2$. Because the conditional probability $P\left( {{n}_{2}}|{{n}_{1}} \right)$ is equal to the fraction of ${{n}_{2}}>{{n}_{1}}$ in the domain $\left\{ \left( {{n}_{1}},{{n}_{2}} \right)|0<{{n}_{1}}<\tau /n,0<{{n}_{2}}<\tau /n \right\}$, ${{n}_{1}}{{n}_{2}}$ take the value ${{h}^{2}}$ with probability $\frac{1}{2}\alpha \frac{\tau }{n}$ and the value 0 with the probability $1-\frac{1}{2}\alpha \frac{\tau }{n}$. According to the central limit theorem, a very large collection of realizations ${{n}_{1}}{{n}_{2}}$ is given by a Gaussian distribution with mean $\frac{1}{2}\alpha {{h}^{2}}\tau$. This result is analogous to $\frac{1}{2}{{\sigma }^{2}}W_{t}^{2}$ which was obtained from the Stratonovich calculus integral of $\mathop{\int }_{0}^{t}{{W}_{t}}\text{d}{{W}_{t}}$. Moreover, the factor of $1/2$ originates in the time order of ${{n}_{1}}$ and ${{n}_{2}}$ which is traced to the constraint ${{n}_{1}}<{{n}_{2}}$ when $P\left( {{n}_{2}}|{{n}_{1}} \right)$ is calculated. In other words, the present value ${{n}_{1}}$ is influenced in a time-symmetric manner and this factor of $1/2$ reflects this symmetry through a properly ordered time integration.

In this simple argument one sees that the δ-autocorrelated noise (used in the definition of white noise) and the Stratonovich calculus approach have the same origin. Both formalisms use different approaches to capture the accumulated influence of short time-scale correlations of a noise source that are not represented in the long time-scale dynamics. Importantly, temporal continuity of the noise is not a necessary condition. In order to demonstrate this, we used discontinuous autocorrelated shot noise to recover both the δ-function definition of white noise and the mid-point selection procedure of Stratonovich calculus. Moreover, one can make the simple case more realistic by treating h as a random variable and the same logic holds.

Finally, we note that despite our pedagogical example $\text{d}x/\text{d}t=xF\left( t \right)$, the results are general. Thus, when we replace x in the right hand side with $b\left( x,t \right)\frac{\partial }{\partial x}b\left( x,t \right)$, we can generalize our argument to any multiplicative noise case. There are of course a myriad of ways to decide which calculus is most appropriate to the problem at hand. Hence, whereas the integral Taylor expansion method enables one to only focus on the nature of the short time-scale processes, the comparison between the integral of autocorrelation and the variance of the short time-scale process focuses the choice on the specific mathematical model or scientific problem. Of particular interest to us is the question of how the stability of the non-autonomous SDEʼs of interest in climate dynamics are influenced by such considerations [20].

According to our result using shot noise, even an infinitesimal noise correlation in a stochastic differential equation can be interpreted using Stratonovich calculus. One might consider our analysis as a generalised version of the Wong and Zakai [1] approach because shot noise here is not a ${{\text{C}}_{1}}$ function. Considering the ubiquity of colored noise in real systems, Itô calculus might be interpreted as an idealized mathematical procedure, only applicable to true white noise processes that are never realized in nature. Hence, we ask whether there are situations in which Itô calculus can be used for colored noise?

The conceptual model we are dealing with in stochastic differential equations is that we have a principle deterministic process whose fate is influenced by short time-scale fluctuations, which are considered as noise. In building a mathematical model it is common to ignore the influence of the short time scale processes on the deterministic dynamics, but there are situations when this may be known to be a poor assumption such as is in the presence of inertial [21] or feedback [22] effects. Indeed, Kupferman and colleagues [21] studied systems with multiplicative colored noise and inertia to find that if the correlation time of the noise is faster (slower) than the relaxation time, this leads to the Itô (Stratonovich) calculus form of the limiting stochastic differential equation. Similarly, Itô calculus is invoked to interpret experiments wherein the time delay of the feedback is much larger than the noise correlation time [22].

These results may also be reinterpreted within the framework of our shot noise formalism. When the inertial or feedback time scales are much shorter than the noise decorrelation time, one can ignore the former and focus on the latter using the logic that lead us to the Stratonovich calculus interpretation. However, when the inertial or feedback time scales are much longer than the noise decorrelation time, the duration w of the shot noise is no longer equivalent to the decorrelation time scale of the noise. Rather, in order for the principle dynamics to be valid, we must interpret w as the inertial or feedback time scale. Hence, during a time increment w, we must assume that two consecutive random events ${{n}_{1}}$ and ${{n}_{2}}$ are independent, $P\left( {{n}_{1}},{{n}_{2}} \right)=P\left( {{n}_{1}} \right)P\left( {{n}_{2}} \right)$, which is analogous to the equation (24). Therefore, despite the colored noise in this case, Itô calculus prevails. Furthermore, this description provides an unambiguous interpretation of the discrete nature of a system [18], which is believed to be a central criterion for the use of Itô calculus.

When a system is described as having multiplicative noise, the criterion for choosing which stochastic calculus is most appropriate—Itô or Stratonovich—has been a source of confusion and great discussion. In most areas of physical science where white noise is defined in terms of a δ-function autocorrelation, Stratonovich calculus is preferred, mainly due to the consistency with the results emerging from the Fokker–Planck equation and the fluctuation-dissipation theorem. In economics, where the Martingale property is considered as the most essential aspect of stochastic random variables, Itô calculus is widely accepted and used in the development of models. On the other hand, in population dynamics, Itô calculus is viewed as a proper continuous approximation of an underlying discrete model in which the Martingale property is guaranteed by construction. However, it would appear more prudent to choose the calculus depending on a clear set of objective considerations.

The core difference between the approaches is seen through the presence of the Martingale property in the multiplicative noise term of the stochastic Langevin equation, and hence manifests itself in the drift term of the associated Fokker–Planck equation. Thus, the principal ambiguity is associated with the ensemble mean of the product of two Brownian motion processes and in the original form it is difficult to assess the validity of the Martingale property. Here we used integral Taylor expansions to pinpoint the source of the deviation between the two calculi, which resides in the characteristics of the short time-scale noise process, the nature of which is described by the autocorrelation. The approach allows one to see the origin of the Stratonovich calculus, which resides in the mid-point selection scheme that is synonymous with including the effect of a finite autocorrelation.

It appears that most realistic signals simulated by white noise have non-Markovian structure, that is, a finite decorrelation time ${{\tau }_{c}}$. When one writes down a model stochastic differential equation, one typically initially assumes that ${{\tau }_{c}}$ is very small relative to the characteristic deterministic dynamical time scales, in which case white noise is a proper approximation. Thus, in such a case, the integral of the noise autocorrelation function and the variance of the noise are of a similar order, and the δ-function autocorrelation is a good approximation that is consistent with Stratonovichʼs calculus. However, when the integral of the noise autocorrelation function is much smaller than the variance of the noise, the δ-function autocorrelation is not an appropriate definition of white noise and Itôʼs calculus is appropriate.

Focusing on the noise itself, in the microscopic limit it is no longer Brownian but is instead represented by discontinuous finite-time signals, or shot noise. A small but finite autocorrelation of shot noise can be characterized by a δ-function, and the accumulation over time of its autocorrelation is represented by the mid-point selection procedure in the Riemann sum of Stratonovich calculus. This demonstrates that the origin of δ-correlated noise and Stratonovich calculus is the infinitesimal autocorrelation of the stochastic noise. Thus, the consistency between the Fokker–Planck equation and Stratonovich calculus has the same physical source, the finite autocorrelation of the stochastic noise that is ignored on time scales long relative to the fluctuation time scales. Therefore, a finite noise autocorrelation is the key criterion for choosing a stochastic calculus, but we must take care in interpreting its origin generally as well as on a case by case basis.

WM thanks NASA for a graduate fellowship and JSW thanks the John Simon Guggenheim Foundation, the Swedish Research Council and a Royal Society Wolfson Research Merit Award for support.