Analyzing a stochastic process driven by Ornstein-Uhlenbeck noise

TL;DR

This paper introduces a direct method to estimate parameters of a non-Markovian stochastic process driven by Ornstein-Uhlenbeck noise, using moments of increments and stochastic Taylor-expansion, simplifying analysis from time series data.

Contribution

It presents a novel, moment-based estimation approach for non-Markovian processes driven by Ornstein-Uhlenbeck noise, avoiding complex Bayesian or embedding techniques.

Findings

Derivation of analytic expressions for moments of process increments.

Establishment of a regression-based parameter estimation method.

Potential for improved analysis of non-Markovian stochastic processes.

Abstract

A scalar Langevin-type process $X(t)$ that is driven by Ornstein-Uhlenbeck noise $\eta(t)$ is non-Markovian. However, the joint dynamics of $X$ and $\eta$ is described by a Markov process in two dimensions. But even though there exists a variety of techniques for the analysis of Markov processes, it is still a challenge to estimate the process parameters solely based on a given time series of $X$. Such a partially observed 2D-process could, e.g., be analyzed in a Bayesian framework using Markov chain Monte Carlo methods. Alternatively, an embedding strategy can be applied, where first the joint dynamic of $X$ and its temporal derivative $\dot X$ is analyzed. Subsequently the results can be used to determine the process parameters of $X$ and $\eta$. In this paper, we propose a more direct approach that is purely based on the moments of the increments of $X$, which can be estimated for different time-increments $\tau$ from a given time series. From a stochastic Taylor-expansion of $X$, analytic expressions for these moments can be derived, which can be used to estimate the process parameters by a regression strategy.

Equations

$$\frac{\partial }{\partial t}X$$

$$\frac{\partial }{\partial t}\eta$$

$$\displaystyle\Delta X(\tau,t_{0})\big{|}_{x_{0}}$$

$$\eta (t)$$

$$\displaystyle\big{<}..\big{>}_{\xi,\eta(0)}$$

$$V(t)$$

$$⟨\eta (t)\eta (t^{\prime })⟩$$

$$dX(t)$$

$$dh[X(t)]$$

$$L_0$$

$$L_1$$

$$dh\{t\}$$

$$h\{t\}$$

$$[L_{0}h]\{s\}$$

$$[L_{1}h]\{s\}$$

$$h\{t\}$$

$$\mathbf{\alpha }$$

$$h[X(t)]$$

$$c_{({\alpha }_1,\dots ,{\alpha }_n)}[X(0)]$$

$$J_{({\alpha }_1,\dots ,{\alpha }_n)}(t)$$

$$d{Z}_j(s)$$

$$J_{(0,\dots ,0)}(t)$$

$$\displaystyle X(t)\big{|}_{x_{0}}$$

$$c_{({\alpha }_1,\dots ,{\alpha }_n)}(x_0)$$

$$c_{(0)}$$

$$c_{(1)}$$

$$\displaystyle\Big{<}F\big{[}\eta(t)\big{]}\big{|}_{x_{0}}\Big{>}$$

$$M^{(1)}(\tau ,x)$$

$$M^{(2)}(\tau ,x)$$

$$\displaystyle\Delta X(\tau)\big{|}_{x}$$

Full text

Markov processes, Stochastic processes

I Introduction

A stochastically forced, first order differential equation provides an appropriate description for the evolution of many physical, chemical or biological systems. For simplicity, we restrict ourselves to the evolution of the scalar quantity $X(t)$ in the following. Additionally, we assume that the coefficient functions in its evolution equation do not explicitely depend on time. Thus, we consider an equation of the form

[TABLE]

where $\eta(t)$ denotes the stochastic force. Such an equation arises not only for the ”obvious” case, where a deterministic system is driven by some external stochastic force, but also for complex dynamical systems consisting of a large number of subsystems. Here, the phenomenon of self-organisation can give rise to a dynamic of so-called order parameters that ”enslave” the dynamic of the microscopic subsystems haken04 , leading to an equation of the above form. However, now the stochastical force $\eta(t)$ can no longer be considered to be external but is an intrinsic part of the system dynamic.

So far, the statistical properties of $\eta(t)$ have not been specified. In practice, this force quite often is treated as Gaussian white noise. Frequently, the central limit theorem can be invoked, which then justifies the assumption of a Gaussian probability density. The assumption of delta-correlated noise, however, is an idealization. Real-world systems usually have some finite correlation time $\theta$. How strong the correlations of $\eta(t)$ affect the statistics of $X(t)$ depends on the ratio of $\theta$ and the characteristic time-scale $T$ of $X(t)$ haenggi95 . For $\theta\ll T$ the force $\eta(t)$ can be approximated by delta-correlated noise, leading to a Markovian description. The probably most famous example is given by Einsteins description of Brownian motion by a Wiener-process einstein05 . Even though the true process is non-Markovian on a microscopic scale, the Markov property can be taken for given for increments larger than some limit timescale. This approach has also successfully been applied to other problems like the description of turbulent velocity increments by a process in scale friedrich97 ; renner01 . Here, the limit timescale is replaced by its spatial analogon, which is denoted as Markov-Einstein coherence length in lueck06 .

Although there are many systems where $\theta$ is sufficiently small and can be neglected, in a variety of systems such an idealization leads to notable differences. For such systems, it is no longer justified to ignore the correlations of $\eta(t)$. However, if we want to account for these correlations, we need a description of the evolution of $\eta(t)$ that goes beyond a ”purely random Gaussian process” risken89 . The most natural and simple generalization of Gaussian white noise is provided by exponentially correlated Gaussian noise, as generated by a stationary Ornstein-Uhlenbeck process. Even if this is not the most general description of colored noise, the assumption that $\eta(t)$ is a stationary process obeying

[TABLE]

makes Eq. (1) applicable to a much larger class of problems. Here, $\xi(t)$ denotes Gaussian white noise with $\left<\xi(t)\right>\!=\!0$ and $\left<\xi(t)\xi(t^{\prime})\right>\!=\!\delta(t\!-\!t^{\prime})$. The characteristic time-scale of $\eta(t)$ is determined by the parameter $\theta$. In the limit $\theta\to 0$, the case of Gaussian white noise is recovered.

Based on this assumption for $\eta(t)$, Eqs. (1) and (2) describe a Markov process in two dimensions. However, the analysis of this process is hampered by the fact that usually only a 1D-series of values of $X(t)$ will be available in practice. In the mathematical community this problem is known as ”partially observed diffusions” dembo86 ; campillo89 . There are approaches to deal with such problems. In a Bayesian framework, e.g., one could use Markov chain Monte Carlo methods for an estimation of $f$, $g$ and $\theta$ (see, e.g., golightly08 ). Alternatively, an embedding approach could be used, where first a series of velocities $\dot{X}(t)$ is calculated from the values of $X(t)$ and subsequently the 2D-system $[X(t),\dot{X}(t)]$ is analyzed. The drift- and diffusion functions of this 2D-system can then be used to determine $f$, $g$ and $\theta$. However, some care has to be taken with this latter approach, because the velocities need to be estimated numerically. This leads to spurious correlations that may affect the results lehle15 .

Here, we propose a more direct approach that is purely based on the moments of the conditional increments

[TABLE]

where $(..)|_{x_{0}}$ denotes conditioning on $X(t_{0})=x_{0}$. In the following, we restrict ourselves to a statistically stationary process $X(t)$, i.e., the moments of $\Delta X$ do not depend on $t_{0}$ and can be estimated from a single time series of $X(t)$ for different values of $\tau$ and $x_{0}$. Our strategy for parameter estimation is based on a stochastic Taylor expansion, which allows us to express the increments $\Delta X$ by an infinite sum that involves multiple integrals with respect to a noise generating process. Since the Ornstein-Uhlenbeck process can be solved analytically, explicit expressions for the correlations of these integrals can be given. Assuming $\tau$ and $\theta$ to be small as compared to the characteristical time-scale of $X(t)$, the moments of $\Delta X$ can be approximated by a finite base of functions $r_{i}(\tau,\theta)$ that are weighted by coefficients $\lambda_{i}(x_{0},\theta)$. These functions are used to fit the moments of $\Delta X$ and thus allow for an estimation of the process parameters. In a first step, the parameter $\theta$ is estimated by a non-linear minimization procedure. Subsequently, by linear fitting, the coefficients $\lambda_{i}$ are estimated, which are then used to determine $f$ and $g$.

It may be noted that in the limit $\theta\to 0$ the functions $r_{i}$ reduce to powers of $\tau$. The above approach then recovers the so-called direct estimation method, which is applicable to processes driven by Gaussian white noise. Due to its simplicity of use, this method has found wide-spread use. For an overview see e.g. friedrich11 .

The standard rules for integro-differential equations apply to the calculations in this paper for $\theta>0$, i.e., for correlated driving noise. Therefore, we adopt the Stratonovich definition of stochastic integrals in the limit $\theta=0$. By this choice, the results will remain valid in unchanged form also in the white-noise limit vankampen81 .

This paper is structured as follows. In Sec. II, we compile some properties of the stochastic force $\eta(t)$. Subsequently, the stochastic Taylor expansion of $X(t)$ is given in Sec. III, which is used in Sec. IV to provide a series representation for the moments of the conditional increments of $X$. The functional form of the series terms is discussed in Sec. V. Subsequently, a series truncation is performed in Sec. VI, which then is used in Sec. VII to formulate a strategy for parameter estimation. To verify the analytical results, a numerical example will finally be given in Sec. VIII.

II Stochastic force

We assume that $\eta(t)$ is a stationary Ornstein-Uhlenbeck process obeying Eq. (2). Ornstein-Uhlenbeck processes are well understood and can be solved analytically. The realization of $\eta(t\geq t_{0})$ can explicitely be expressed in terms of an initial value $\eta(t_{0})$ and the realization of $\xi(t\geq t_{0})$. Since we focus in the following on the stationary process, we are free to choose $t_{0}\equiv 0$, which simplifies notation. One then finds

[TABLE]

This equation holds for arbitrary values of $\eta(0)$. It describes a realization of the process $\eta(t)$, i.e., a trajectory in time, in terms of $\eta(0)$ and a trajectory of $\xi(t)$. Expectation values of functionals of $\eta(t)$ thus are obtained by averaging over the realizations of $\eta(0)$ and $\xi(t)$. As these quantities are statistically independent, averaging may be performed in two steps using

[TABLE]

In the next section, $\eta(t)$ will be expressed as derivative of a noise-generating process $V(t)$ with $V(0)\equiv 0$. This implies

[TABLE]

This process plays a comparable role for $\eta(t)$ as the Wiener process does for Gaussian white noise. Using Eq. (4), we also may describe a trajectory of $V(t)$ directly in terms of $\eta(0)$ and a trajectory of $\xi(t)$,

[TABLE]

It may easily be checked that $V(t)$ approaches the Wiener process $W(t):=\int_{0}^{t}\xi(s)\,ds$ in the limit $\theta\to 0$.

Finally, we consider some properties of the stationary process. With Eq. (4) the stationary probability density function (PDF) of $\eta(t)$ is found to be a Gaussian with variance $1/(2\theta)$ and vanishing mean. Furthermore, the autocorrelation function is found to be

[TABLE]

This means that $\eta(t)$ is normalized in the following sense: its strength, i.e., the integral over its autocorrelation function, is constant and equals unity. Therefore, in the limit $\theta\!\to\!0$ the autocorrelation approaches $\delta(t-t^{\prime})$ and $\eta(t)$ approaches Gaussian white noise $\xi(t)$.

III Stochastic Taylor expansion

Since, later on, we will focus on moments of conditional increments of $X(t)$, we need an analytic description of these increments. Assuming smooth functions $f$ and $g$, such a description can be provided by a stochastic Taylor expansion, which allows us to express a trajectory of $X(t)$ in terms of $X(0)$, values and derivatives of $f$ and $g$ at $X(0)$, and a trajectory of $\eta(t)$. Such expansions are described in great detail, e.g., in platen99 , and we will closely follow these lines.

The starting point is Eq. (1) in the form $dX=fdt+g\eta dt$. Expressing $\eta(t)$ as derivative of a noise generating process $V(t)$ as defined by Eq. (6), this may be written as

[TABLE]

Next, the infinitesimal increment of an arbitrary, smooth, function $h(X)$ is considered. Since $X(t)$ is continuously differentiable, the standard chain-rule of differentiation applies,

[TABLE]

Expressing $dX$ by Eq. (10) and introducing the operators

[TABLE]

this may be written as

[TABLE]

Here, we use the notation $\{t\}$ to indicate that all arguments of a function or expression are to be evaluated at time $t$. Now the actual expansion of $h$ can be started. In integral form, Eq. (13) reads

[TABLE]

Since we considered $f$, $g$ and $h$ to be smooth, $L_{0}h$ and $L_{1}h$ are also smooth functions of $X$. Consequently, Eq. (15) can be applied,

[TABLE]

Inserting these results into Eq. (15) then yields

[TABLE]

The single integrals from Eq. (15), which had time-dependent integrands $[L_{i}h]\{s\}$, are now replaced by single integrals with constant integrands $[L_{i}h]\{0\}$ plus additional double integrals with time-dependent integrands $[L_{i}L_{j}h]\{s^{\prime}\}$. Expressing these functions by Eq. (15) will put the game on the next level, leading to constant double integrals plus variable triple integrals — and so on. In the end, one is left with an infinite sum of multiple integrals, which only depend on $t$ and the realization of $V(t)$, that are multiplied by coefficient functions that only depend on values and derivatives of $f$, $g$ and $h$ at $X(0)$,

[TABLE]

Using a multi-index ${\boldsymbol{\alpha}}$, defined as

[TABLE]

the expansion of $h[X(t)]$ can be compactly written as

[TABLE]

where the coefficient functions $c_{\boldsymbol{\alpha}}$ are given by

[TABLE]

and the integrals $J_{\boldsymbol{\alpha}}$ by

[TABLE]

with

[TABLE]

In general, these integrals are functionals of the realization of $V(t)$ respectively $\eta(t)$ and thus stochastic quantities. Only for $\alpha_{1}=\ldots=\alpha_{n}=0$ the integrals become purely deterministic and evaluate to

[TABLE]

So far, the expansion of some arbitrary function $h(X)$ has been considered. Being interested in the expansion of $X(t)$ itself, we choose $h(X)\equiv X$ in the following. Additionally, we fix the value $X(0)$ to $x_{0}$, which then leaves us with

[TABLE]

where the coefficient functions are now defined as

[TABLE]

Omitting arguments and using a prime to denote derivatives with respect to $X$, the first few of these functions (to be evaluated at $x_{0}$) read

[TABLE]

The conditioning of $J_{\boldsymbol{\alpha}}$ in Eq. (39) deserves some comment. After all, a realization of $J_{\boldsymbol{\alpha}}(t)$ does not depend on $X(0)$ but is a pure functional of $\eta(t)$, which itself is a functional of $\eta(0)$ and the realization of $\xi(t)$. However, the PDF of $\eta(0)$ will, in general, depend on $X(0)$ [see, e.g., Appendix D for the expectation value of $\eta(0)|_{x_{0}}$]. As a consequence, ensemble averages of any conditioned functional $F[\eta(t)]|_{x_{0}}$ need to be calculated by averaging over the realizations of $\xi(t)$ and over the conditional realizations $\eta(0)|_{x_{0}}$. Averaging may still be performed in two steps, e.g., by

[TABLE]

IV Conditional moments of $\Delta X$

We now turn to mean and variance of the conditional process increments of $X(t)$,

[TABLE]

where the increments are denoted by

[TABLE]

Since we are conditioning on the value of $X$ at some arbitrary time $t$, we denote this value by $x$ instead of by $x_{0}$. Additionally, we suppress the function argument $t$, because the statistical properties of the increments $\Delta X$ do not depend on time for a stationary process. Stationarity also implies that the moments $M^{(k)}$ can be estimated from a given time series of $X(t)$ by replacing the above ensemble-averages by time-averages (tacitly assuming ergodicity). Using the results from the previous section, we already have an analytical description for the increments,

[TABLE]

Hence, the conditional moments are given by

[TABLE]

with (omitting arguments)

[TABLE]

As a result, we now have analytic descriptions of the moments — but unfortunately in terms of infinite series. In order to obtain approximate descriptions with a finite number of terms, the functional form of $\phi_{\boldsymbol{\alpha}}(\tau,x)$ needs to be investigated in the following. This also provides us with the functional form of $\phi_{{\boldsymbol{\alpha}},{\boldsymbol{\beta}}}(\tau,x)$, because a product $J_{\boldsymbol{\alpha}}J_{\boldsymbol{\beta}}$ can be expressed by a sum of integrals $J_{\boldsymbol{\gamma}}$ (see Appendix A),

[TABLE]

which implies

[TABLE]

V Functional form of $\phi_{\boldsymbol{\alpha}}$

The starting point for the calculation of $\phi_{\boldsymbol{\alpha}}$ is the definition of the integral $J_{\boldsymbol{\alpha}}$, as provided by Eq. (34). Let us consider an index vector ${\boldsymbol{\alpha}}$ of length $n$ and denote the number of its non-zero entries by $m$. Using $dV=\eta\,dt$, Eq. (34) may then be written as $n$-fold integral with respect to time over an $m$-fold product of $\eta$,

[TABLE]

Here, the shortcut $\Omega(\tau)$ has been introduced to denote the integration domain (a simplex in ${\mathbb{R}}^{n}$ with $0\leq s_{i}\leq s_{i+1}$ and $s_{n}\leq\tau$). The ensemble average of J_{\boldsymbol{\alpha}}\big{|}_{x} then reads

[TABLE]

where the values $j_{1},\ldots,j_{m}$ denote the positions of the non-zero entries in ${\boldsymbol{\alpha}}$ and $C_{\eta}$ the $m$-point correlation function of \eta(t)\big{|}_{x}. For arbitrary times $t_{1},\ldots,t_{m}$ this function is defined as

[TABLE]

According to Eq. (4), $\eta$ may be splitted up into one part depending only on $\eta(0)$ and another one depending only on $\xi$,

[TABLE]

with

[TABLE]

Consequently, $C_{\eta}$ can be expressed in terms of conditional moments of $Y$ and correlation functions of $u(t)$, denoted as

[TABLE]

For example, we find

[TABLE]

Explicit expressions for $C$ can be found by virtue of Eq. (56). It turns out that the correlation functions of $u(t)$ have the same structure as those of Gaussian white noise $\xi(t)$ (see Appendix B). The $k$-point correlation of $u$ vanishes for odd values of $k$, while for even values it can be expressed by a sum of products of the two-point correlation [Eq. (133)].

Since $C(t_{1},\ldots,t_{k})$ vanishes for odd $k$, the expressions for $C_{\eta}$ may contain either only even or only odd moments of Y\big{|}_{x}. To provide an example:

[TABLE]

With the above results, the intgral on the right-hand side of Eq. (52) can be evaluated, which leads to

[TABLE]

Here, we introduced the shortcuts $a_{k}(\tau)$ to denote the functions that stem from the integrations with respect to time. Actually, these functions depend also on the index vector ${\boldsymbol{\alpha}}$, but we supressed this argument for notational simplicity. For an explicit example see Appendix C.

Later on, it proves to be useful to re-arrange the right-hand side of this equation by expressing the powers of $Y$ in terms of Hermite polynomials in $Y$. Reordering terms then yields

[TABLE]

where the $k$-th Hermite polynomial is defined as

[TABLE]

and the functions $b_{k}$ are linear combinations of the functions $a_{k}$. For example, a right-hand side of the form \big{<}Y^{2}|x\big{>}a_{0}+a_{1} becomes \big{<}H_{2}(Y)|x\big{>}b_{0}+b_{1} with $b_{0}=a_{0}/4$ and $b_{1}=a_{1}+a_{0}/2$.

By mathematical induction, it may be shown that the functions $a_{k}(\tau)$, and thus also the functions $b_{k}(\tau)$, are linear combinations of the functions

[TABLE]

with

[TABLE]

and

[TABLE]

We thus may express $b_{k}$ in the form $\sum\tilde{\lambda}_{ij}\tilde{r}_{ij}$. Note, that the coefficients $\tilde{\lambda}_{ij}$ do not depend on $\theta$, because this dependency is completely accounted for by the functions $\tilde{r}_{ij}$. This property will be useful in the next section, when we consider the magnitude of individual terms. However, since the functions $\tilde{r}_{ij}$ depend on $\ell({\boldsymbol{\alpha}})$, this property can not be sustained for the following base of functions $r_{ij}$, which is used to describe the $\tau$-dependency of $b_{k}$, and thus of $\phi_{\boldsymbol{\alpha}}$, for arbitrary vectors ${\boldsymbol{\alpha}}$,

[TABLE]

with

[TABLE]

Note that the product of any two functions of this base lies in the linear span of $\cal B$. According to Eq. (50), therefore, $\cal B$ not only provides a base for the functions $\phi_{\boldsymbol{\alpha}}$ but also for the functions $\phi_{{\boldsymbol{\alpha}}{\boldsymbol{\beta}}}$.

VI Series truncation

With the results from the previous section, the series representation of the moments $M^{(k)}$, Eq. (47), can now be expressed in terms of functions $r_{i}\in{\cal B}$,

[TABLE]

where each coefficient function $\lambda_{i}$ consists of an infinite sum of terms. These terms, in general, are formed by powers of $\theta$, the functions $c_{\boldsymbol{\alpha}}$ from the Taylor expansion [Eq. (40)], and the expectation values $\left<H_{n}(Y)|x\right>$.

In order to approximate $M^{(k)}$ by a finite number of functions, it becomes necessary to make some assumptions on the magnitude of the individual terms. First, we assume that $X(t)$ has been normalized to ensure a characteristic time-scale of unity and coefficient functions $c_{\alpha}$ of order $O(1)$. Second, we assume

[TABLE]

where $\varepsilon$ has been introduced to denote a quantity that is small as compared to unity. We may identify $\varepsilon$ as the largest increment $\tau$ for which we assume that our truncated description of $M^{(k)}$ holds. The terms $\left<H_{n}(Y)|x\right>$, finally, are (for the time being) treated as $O(1)$ terms, because $Y:=\theta^{1/2}\eta$ is a Gaussian random variable with a constant variance of $1/2$. It remains to ask for the magnitude of $r_{i}(\tau)$. According to Eq. (71), this is a term of order $O(\varepsilon^{a})$ for $r_{i}=r_{a0}$, whereas for $r_{i}=r_{0b}$ or $r_{i}=r_{ab}$ it may be treated as term of order $O(1)$, because in this case the value range of $r_{i}$ is finite and depends only on $a$ and $b$.

We now focus on a description of $M^{(1)}$, in which only terms up to order $O(\varepsilon^{3})$ are considered. According to Eq. (66) and the above assumptions, a function $\phi_{\boldsymbol{\alpha}}$ may only give rise to terms of order $O(\theta^{j}\tau^{k})$ with $j+k=\ell({\boldsymbol{\alpha}})$. Therefore, the lowest order contributions in terms of $\varepsilon$ are of order $O(\varepsilon^{\ell({\boldsymbol{\alpha}})})$. We thus may write Eq. (47a) in the form

[TABLE]

Evaluating these functions $\phi_{\boldsymbol{\alpha}}$ and re-sorting terms then provides an approximation of $M^{(1)}$ in terms of seven base functions $r_{i}$. Of course, the resulting coefficients of these functions are only truncated versions of the coefficients $\lambda_{i}$ in Eq. (72). The coefficient of $r_{1,0}\equiv\tau$, e.g., is only accurate up to order $O(\theta)$ — but this will be sufficient for our purposes.

So far, the expectation values $\left<H_{n}(Y)|x\right>$ have been treated as terms of order $O(1)$. Actually, however, this is only a lower limit for their order of magnitude. This becomes obvious by looking at the Fokker-Planck equation of the stationary 2D-process $[X(t),\eta(t)]$, where it turns out that $\left<H_{1}(Y)|x\right>$ is of order $O(\varepsilon)$ (see Appendix D),

[TABLE]

For $n>1$, such explicit results are not available. Nevertheless, the magnitude of terms can be shown to obey (see Appendix E)

[TABLE]

With these findings, a number of terms become sufficiently small to be neglected, which leads to an approximation of $M^{(1)}$ in terms of the function base $\{r_{0,1},r_{1,0},r_{2,0},r_{3,0}\}$. Additionally, it turns out that terms $\left<H_{n}(Y)|x\right>$ with $n>1$ are no longer present in the coefficients of these functions.

As a last step, we switch to a modified base $\{r_{0,1},r_{1},r_{2},r_{3}\}$, where the new base functions $r_{i}$ are linear combinations of the former ones,

[TABLE]

This base not only leads to simpler coefficients in general, but most importantly, the coefficient of $r_{0,1}$ now becomes sufficiently small to be neglected, which leaves us with a base of only three functions.

Following the above lines, also an approximation of $M^{(2)}$ can be obtained. As it is the case for the approximation of $M^{(1)}$, calculations are straightforward but cumbersome. Therefore, we only give the final results here, which can be summarized as follows. The moments $M^{(k)}$ can be approximated by

[TABLE]

with

[TABLE]

The coefficients of $r_{1}$ are found to be (omitting arguments)

[TABLE]

These equation will allow us to determine $f$ and $g$, once we manage to provide values for $\lambda^{(k)}_{1}$ and $\theta$.

VII Parameter estimation

Now the estimation of $\lambda^{(k)}_{i}$ and $\theta$ can be addressed. For a given time series of $X$, the moments $M^{(k)}$ can be estimated for a number of $N$ time-increments $\tau_{\nu}$ with

[TABLE]

These estimates of $M^{(k)}$ can be fitted by means of Eq. (80) in a least-square sense. Estimates of $\theta$ and $\lambda^{(k)}_{i}$ thus may be obtained by minimizing the residuals

[TABLE]

where the values $\lambda^{(k)}_{i}$ have been combined into the vector ${\boldsymbol{\lambda}}^{(k)}$ for syntactical convenience. Additionally, the dependency of the functions $r_{i}$ on $\theta$ has been made explicit by the syntax. Due to this dependency, a non-linear approach is needed for the minimization of the above residuals.

It may be noted that minimizing $R^{(k)}$ includes the estimation of $\theta$ for each value of $x$ and $k$ — despite the fact that $\theta$ is a constant. This is neither the most efficient nor the most accurate way for parameter estimation. Instead, we will estimate $\theta$ only once, based on the autocovariance $A(\tau):=\left<X(\tau)X(0)\right>$, respectively its increments

[TABLE]

Estimates of $\Delta A$ are much more accurate than estimates of $M^{(k)}$, because the latter are based on far less data, due to the conditioning on $x$. An approximation of $\Delta A$ that is accurate up to terms of order $O(\varepsilon^{3})$ is given by (see Appendix G)

[TABLE]

An estimate of $\theta$ may thus be obtained by minimizing

[TABLE]

For fixed $\theta$, the optimal values $\lambda_{i}^{*}(\theta)$ can be obtained by linear regression. This means: we can explicitly calculate

[TABLE]

as well as the corresponding residual value $R[{\boldsymbol{\lambda}}^{*}(\theta),\theta]$. The optimal value $\theta^{*}$, which corresponds to the global minimum $R[{\boldsymbol{\lambda}}^{*}(\theta^{*}),\theta^{*}]$, is thus formally given by

[TABLE]

In practice, $\theta^{*}$ may be found numerically, e.g., by using some recursive strategy to search for the minimum of $R[{\boldsymbol{\lambda}}^{*}(\theta),\theta]$ within the interval $[0,\theta_{\text{max}}]$. We safely may choose $\theta_{\text{max}}=\tau_{\text{max}}$, because, according to our assumptions on the magnitude of $\tau$ and $\theta$, Eq. (73), we anyway need to rely on $\theta<\tau_{\text{max}}$. Otherwise our series-truncation, Eq. (80), would no longer be valid.

Once $\theta$ has been estimated, estimates of $\lambda^{(k)}_{i}$ become accessible by a linear regression strategy, which allows us to explicitly calculate

[TABLE]

As a last step, it remains to determine $f$ and $g$ using the estimates of $\theta$ and $\lambda^{(k)}_{1}$. This is achieved by writing Eq. (84) in the form (omitting arguments and dropping asteriscs)

[TABLE]

Because $\theta$ is assumed to be small, $f$ and $g$ can be determined by a fixed-point iteration. For given values $f^{(n)}$ and $g^{(n)}$, the right-hand sides of the above equations provide the definitions for $f^{(n+1)}$ and $g^{(n+1)}$. However, as this requires the evaluation of spatial derivatives of $f^{(n)}$ and $g^{(n)}$, one needs to simultaneously iterate the values at different locations $x_{i}$. The required derivatives can then be estimated by some numerical differencing scheme. Appropriate starting values are provided by $f^{(0)}(x_{i})=\lambda^{(1)}_{1}(x_{i})$ and $g^{(0)}(x_{i})=[\lambda^{(2)}_{1}(x_{i})]^{1/2}$ .

VIII Numerical example

In the following, we investigate a numerical test case, for which we use a non-symmetric, heavy-tailed, process with multiplicative noise,

[TABLE]

with

[TABLE]

We use this system of equations to generate discrete time series of $X(t)$, consisting of $10^{7}$ points, using a sampling timestep $dt\!=\!0.005$. Integration is performed using the Euler-scheme with an internal timestep $\delta t\!=\!0.02\times\min(\theta,dt)$. The global time-scale of $X(t)$ can be estimated from its autocorrelation function and approximately equals unity for small values of $\theta$. In Fig. 1 excerpts of the generated time series are shown for different values of $\theta$, and in Fig. 2 the corresponding probability densities $p(X)$ and the increments $\Delta A$ of the autocorrelation functions are provided. While stronger correlations of the driving noise lead to notably smoother time series, almost no effect on the probability density can be seen.

For the analysis of a given time series, the regression functions $r_{i}(\tau,\theta)$ play a central role. Therefore, we give them explicitely here again,

[TABLE]

The actual analysis can be summarized as:

Estimate the correlation time $\theta$ by non-linear fitting the increments of the autocovariance of $X$ with the functions $r_{i}(\tau,\theta)$.

• 2)

Use the estimated value $\theta^{*}$ to estimate the values $\lambda_{1}^{(k)}(x,\theta^{*})$ by linear fitting the moments $M^{(k)}(\tau,x)$ with the functions $r_{i}(\tau,\theta^{*})$.

• 3)

Use the estimated values $\lambda_{1}^{*(k)}(x,\theta^{*})$ to calculate estimates for $f$ and $g$ using Eq. (94).

These steps will now be detailed. We first consider the estimation of the correlation time $\theta$. As mentioned in Sec. VII, an estimate $\theta^{*}$ may be found by minimizing the residual $R[{\boldsymbol{\lambda}}^{*}(\theta),\theta]$, where the vector ${\boldsymbol{\lambda}}^{*}(\theta)$ is obtained from a linear fit of $\Delta A(\tau)$ using the functions $r_{i}(\tau,\theta)$. To find the minimum of $R$ in an interval $[\theta_{\text{min}},\theta_{\text{max}}]$, we use a recursive strategy. First, the residual is evaluated for a number of equidistant values $\theta_{i}$ covering the whole interval. Next, the interval is narrowed and repositioned such that it only covers the vicinity of the value $\theta_{i}^{*}$, for which the residual was found to be smallest. These steps can now be repeated until the desired numerical accuracy is reached.

In our example, we first use the values $\Delta A(\nu\,dt)$ with $1\leq\nu\leq 60$ for the fits. This corresponds to a maximum time increment $\tau_{\text{max}}=0.3$. Since we use a truncated series representation for the description of $\Delta A$, the value of $\tau_{\text{max}}$ affects the systematic errors of the fits and should be choosen as small as possible. Therefore, once we have calculated $\theta^{*}$ with $\tau_{\text{max}}=0.3$, we restrict the maximum increment to $\tau^{*}_{\text{max}}=\sqrt{\theta^{*}}$, which is consistent with our assumptions on the magnitude of terms, and repeat the calculation of $\theta^{*}$.

Estimates for $\theta$ that are obtained by following this strategy are shown in Fig. 3 for the range $0.001\leq\theta\leq 0.1$. Even if $\theta$ seems to be slightly underestimated for $\theta>0.01$, the overall accuracy is quite good.

With an estimate $\theta^{*}$ at hand, the coefficients $\lambda_{1}^{(k)}$ are obtained from linear fits of the moments $M^{(k)}$ using the functions $r_{i}(\tau,\theta^{*})$. For the estimation of these moments, we use a binning approach, where the range $-2\leq x\leq 3$ is divided into 25 bins. For each bin we estimate the values $M^{(k)}(\nu\,dt,x)$ with $1\leq\nu\leq 60$ from the data. Here, $x$ is taken to be the position of the bin-center. In Fig. 4, estimated values and resulting fits of $M^{(k)}$ at $x=-0.9$ are shown for different values of $\theta$. The values obtained from the data can excellently be fitted with the functions $r_{i}$. The mean error is only about $2.5\times 10^{-4}$.

Finally, we consider the estimates of the coefficients $\lambda_{1}^{(k)}$ and of the functions $f$, $f+gg^{\prime}/2$ and $g$. We do not show the results for $\theta=0.001$, because these would look almost identical to the results for $\theta=0.01$, which are presented in Fig. 5. Here, we find that the estimates of $f$ and $g$ are in very good accordance with the true values. Additionally, it shows that — for the given value of $\theta$ — the values of $\lambda_{1}^{(1)}$ and $[\lambda_{1}^{(2)}]^{1/2}$ are almost identical to the values of $f+gg^{\prime}/2$ and g, respectively. But this will change, when larger values of $\theta$ are considered.

In Fig. 6 the results for $\theta=0.05$ are shown, where $[\lambda_{1}^{(2)}]^{1/2}$ clearly deviates from $g$. However, since we account for this deviation by Eq. (94), we still obtain accurate estimates for $f$ and $g$.

To also show the limitations of our approach, the results for $\theta=0.1$ are presented in Fig. 7. Here, the estimates become less accurate outside the range $-1<x<2$. Especially for $x>2$ the estimates of $\lambda_{1}^{(1)}$ now show fluctuations that hamper the estimation of the spatial derivatives that are needed to calculate $f$ and $g$ using Eq. (94).

IX Conclusions

A parameter-free approch has been developed that allows for the analysis of a stochastic process $X(t)$ that is driven by exponentially correlated, Gaussian noise. This analysis is purely based on the moments of the conditional increments of $X$ and provides estimates for the drift- and diffusion-functions of the process as well as for the correlation time $\theta$ of the driving noise.

It should be noted that we use a perturbative approach, where $\theta$ is assumed to be small as compared to the characteristic time-scale of $X$. Actually, the method presented in this paper is accurate up to first order terms in $\theta$. In principle, however, also higher order approximations are possible. The method may be seen as generalization of the direct estimation method friedrich11 , which also formally is recovered in the limit $\theta\to 0$.

The applicability and accuracy of our approach has been demonstrated by a numerical example, where reasonable accurate results are obtained even for values of $\theta$ as large as ten percent of the global time-scale. For smaller values of $\theta$, the results are (aside finite-size fluctuations) close to exact.

The presented approach is straightforward to implement and neither demanding with regard to memory nor to CPU power. An analysis of a series of $10^{7}$ values is performed within a few seconds on a standard desktop PC.

It would be interesting to also apply our approach to an analysis of turbulent velocity increments as described in friedrich97 ; renner01 , where the description of the increments is approximated by a Markov-process in scale. The observed conditional moments of the process-increments show a striking similarity to the moments $M^{(k)}$ [Fig. 4] of a process that is driven by Ornstein-Uhlenbeck noise. However, since the drift- and diffusion functions for this problem clearly depend on scale, our approach first needs to be extended to non-stationary processes, which is a task for the future.

Appendix A Products of integrals $J_{\boldsymbol{\alpha}}$

According to Eq. (34), the integrals $J_{\boldsymbol{\alpha}}$ are defined as

[TABLE]

with

[TABLE]

There are a number of obvious solutions, like

[TABLE]

Additionally, it is possible to express a product of two integrals by a sum of single integrals: According to the above definition, an integral $J_{\boldsymbol{\alpha}}$ may be written as

[TABLE]

where the syntax $(\alpha_{1},\ldots,\alpha_{n})^{-}\!:=\!(\alpha_{1},\ldots,\alpha_{n-1})$ has been introduced. The differential increment of $J_{\boldsymbol{\alpha}}$ is thus given by

[TABLE]

For the differential increment of a product $J_{\boldsymbol{\alpha}}J_{\boldsymbol{\beta}}$ one finds

[TABLE]

In integrated form, we obtain (assuming ${\boldsymbol{\beta}}$ to have $m$ components)

[TABLE]

This equation can now be applied recursively to the products within the integrals. In the end, this leads to a sum of integrals $J_{\boldsymbol{\gamma}}$ with index-vectors of length $n+m$,

[TABLE]

Actually, the vectors ${\boldsymbol{\gamma}}\in{\cal M}({\boldsymbol{\alpha}},{\boldsymbol{\beta}})$ represent all ${n+m\choose n}$ possibilities to mix the indices of ${\boldsymbol{\alpha}}$ and ${\boldsymbol{\beta}}$ while keeping their relative ordering, This means, the position of $\alpha_{i}$ in ${\boldsymbol{\gamma}}$ must always preceed that of $\alpha_{i+1}$, and the same must hold for the components of ${\boldsymbol{\beta}}$. As an example, one obtains

[TABLE]

With this rule at hand, one finds, e.g.,

[TABLE]

which implies (assuming a vector of length $n$)

[TABLE]

One may also derive relations like

[TABLE]

Appendix B Correlation functions of $u(t)$

Using the definition of $u(t)$ [Eq. (56)], the definition of $C$ [Eq. (57)] reads

[TABLE]

As $\xi(t)$ is Gaussian white noise, the expectation values $\left<\xi(t_{1})\cdots\xi(t_{n})\right>$ are well known. For odd values of $n$ they are vanishing,

[TABLE]

All even correlations can be expressed in terms of two-point correlations. For $n=2k$, this leads to a sum of $k$-fold productes of delta functions. This sum contains $1\times 3\times\cdots\times(n\!-\!1)$ terms, which is the number of possibilities to permutate the function arguments of such a product when only distinguishable functions are allowed (functions may be indistinguishable due to the symmetry of the delta function or due to the commutativity of multiplication). Up to $n=4$ this reads

[TABLE]

According to Eq. (122), it follows immediately that $C$ also vanishes for odd values of $n$,

[TABLE]

For $n=2$ one finds

[TABLE]

The higher order correlation functions of $u(t)$ can be expressed in terms of two-point correlations like in the case of Gaussian white noise. This can be seen when inserting the expressions for $\left<\xi(t_{1})\cdots\xi(t_{2k})\right>$ into Eq. (122). For each of the $k$-fold products of the delta function, the integral factorizes into a $k$-fold product of integrals of the form of Eq. (133). Therefore, the structure of the correlation functions of $\xi(t)$ directly translates to that of $u(t)$. One finds

[TABLE]

Appendix C Explicit example for $\phi_{\boldsymbol{\alpha}}$

To provide an explicit example for the calculation of $\phi_{\boldsymbol{\alpha}}$, the case ${\boldsymbol{\alpha}}=(1,0,1)$ is considered. Equation (52) then reads

[TABLE]

where we again denote the integration domain $0\leq s_{1}\leq s_{2}\leq s_{3}\leq\tau$ by $\Omega(\tau)$ and $C_{\eta}$ is defined according to Eq. (53).

By expressing $\eta$ in terms of $Y$ and $u$ [Eqs. (54), (55) and (56)], the two-point correlation $C_{\eta}(s_{1},s_{3},x)$ can be expressed in terms of moments of $Y$ and correlation functions uf $u$, which are denoted by $C$ and defined according to Eq. (57)

[TABLE]

Taking into account that, according to Eq. (125), the functions $C(s_{1})$ and $C(s_{3})$ are vanishing, this leaves us with (omitting index vector and arguments for $\phi$)

[TABLE]

Inserting Eq. (133) and taking into account that $s_{1}\leq s_{3}$ holds within the integration domain, this becomes

[TABLE]

which evaluates to

[TABLE]

with

[TABLE]

In terms of the Hermite Polynomials $H_{0}(Y):=1$ and $H_{2}(Y):=4Y^{2}-2$ this can be re-written as

[TABLE]

with

[TABLE]

Appendix D Expectation value of \eta\big{|}_{x}

Equations (1) and (2) describe a Markov process in two dimensions. Using $x$ and $s$ to denote the phase-space variables of $X$ and $\eta$, the Kramers–Moyal coefficients of the corresponding Fokker–Planck equation are given by

[TABLE]

The Fokker–Planck equation thus reads

[TABLE]

Integrating with respect to $s$ then gives [using $p(x;s)\!=\!p(x)p(s|x)$ and $\int_{s}\!s\,p(s|x)\!=\!\left<\eta|x\right>$]

[TABLE]

where $j$ denotes the probability flux. For the stationary process we have $\partial_{t}p(x)=0$, implying a constant flux. For natural boundary conditions (vanishing flux and density at $|x|\to\infty$) this implies $j\equiv 0$ and thus

[TABLE]

Appendix E Expectation values \big{<}H_{n}(Y)|x\big{>}

In terms of $Y(t)\equiv\sqrt{\theta}\eta(t)$ our evolution equations [Eqs. (1) and (2)] read

[TABLE]

The corresponding Fokker-Planck equation for the stationary process may then be written as

[TABLE]

where $p(x;y)$ denotes the stationary joint PDF of $X(t)$ and $Y(t)$. Actually, this density also depends on the parameter $\theta$, but we will not make this explicit by the syntax. In the following, we express $p(x;y)$ by a Hermite expansion of the form

[TABLE]

with

[TABLE]

The function $p_{0}(x)$, finally, denotes the density of $X(t)$ in the limit $\theta\to 0$. In this limit, Eq. (1) becomes $\dot{X}=f+g\xi$ (to be interpreted in the Stratonovich sense). The stationary density of $X$ can then be calculated from the corresponding Fokker-Planck equation $0=\partial_{x}[fp_{0}-\frac{1}{2}g\partial_{x}(gp_{0})]$, leading to

[TABLE]

where $N$ is a normalization constant.

As the Hermite polynomials $H_{n}(y)$ are orthogonal under the weight $G(y)$, i.e.,

[TABLE]

one first finds from Eq. (155)

[TABLE]

Using $p(x;y)=p(x)p(y|x)$ and $\int_{y}F(y)p(y|x)=\left<F(Y)|x\right>$, the left-hand side of this equation may be rewritten to obtain

[TABLE]

For the case $n=0$, one finds (because of $H_{0}(y)\equiv 1$)

[TABLE]

Together with Eq. (156) and $p_{0}(x)=\lim_{\theta\to 0}p(x)$, this provides us with the value of $c_{0,0}$,

[TABLE]

For $n=1$, a relation between the coefficients $c_{0,i}$ and $c_{1,i}$ can be obtained by using $\left<H_{1}(Y)|x\right>=-2\sqrt{\theta}f/g$ [Eq. (75)] and $p=p_{0}c_{0}$. One first finds

[TABLE]

and as this equation holds for arbitrary $\theta$, Eq. (156) implies

[TABLE]

In the general case one obtains

[TABLE]

Assuming the coefficients $c_{i,j}$ to be of order $O(1)$ then yields

[TABLE]

which implies $\left<H_{n}(Y)|x\right>=O(\theta^{n/2})$, as claimed by Eq. (76).

It remains to be shown, however, that there exists a set of finite coefficients $c_{i,j}$, for which Eq. (155) is a solution of the Fokker-Planck equation as specified by Eq. (154). In order to calculate these coefficients, we insert Eq. (155) into Eq. (154). Using the well known relations

[TABLE]

and

[TABLE]

this first leads to (omitting arguments)

[TABLE]

Multiplying by $H_{n}(y)$ and integrating with respect to $y$ then gives (dividing by $\theta^{n/2}$ and formally defining $c_{-1}:=0$)

[TABLE]

As Eq. (159) implies $\partial_{x}(gp_{0})=\frac{2f}{g}p_{0}$, we can get rid of the factor $p_{0}$. In terms of the operator

[TABLE]

this leads to

[TABLE]

For $n=0$, this equation does not provide any additional information, because it evaluates to $0=L[(f/g)c_{0}+c_{1}]$, which, according to Eq. (166), is fulfilled for all values of $\theta$. Therefore, we only need to look at $n>0$ in the following. Inserting Eq. (156) and sorting terms by powers of $\theta$ then yields

[TABLE]

As this equation must hold for arbitrary $\theta$, all expressions in curly brackets must vanish individually. The first of these expressions, together with $c_{0,0}\equiv 1$, allows us to calculate all coefficients $c_{n,0}$,

[TABLE]

Explicitely one finds

[TABLE]

Similarly, all coefficients $c_{n,1}$ (and subsequently $c_{n,2}$, $c_{n,3}$, $\ldots$) can be calculated using

[TABLE]

However, to start the iterative calculation of $c_{1,i}$, $c_{2,i}$, $\ldots$, we need the coefficient $c_{0,i}$, which can be obtained as follows. For $n=1$, we may use Eq. (167) to epress the left-hand side of the above equation by $-(f/g)c_{0,i}$. This leads to

[TABLE]

We thus find

[TABLE]

with

[TABLE]

where $C_{i}$ is an integration constant, which can be determined by using the fact that $\int_{x}p_{0}c_{0,i}$ vanishes for $i>0$ (see Appendix F). Multiplying Eq. (186) by $p_{0}$ and integrating with respect to $x$, therefore, yields

[TABLE]

To summarize results: We now have equations for all coefficients $c_{i,j}$ and for all integration constants $C_{i}$. But, as noted above, these quantities need to be finite to ensure the validity of Eq. (76). We thus need to presume smooth and finite functions $g$ and $f/g$. Additionally, the limit density $p_{0}$ needs to decay sufficiently fast, to ensure finite values $C_{i}$.

As a final remark: The result for the coefficient $c_{0,1}$, which is found to be

[TABLE]

may be checked for correctness using one of the small-$\theta$ approximations for $p(x)$ that are available in the literature (see, e.g., haenggi95 ). These approximations are known to correctly account for the first order terms in $\theta$. Therefore, when expanding one of them into a power-series in $\theta$, the first order term needs to equal $\theta p_{0}c_{0,1}$ — which indeed is found to be the case.

Appendix F Integrals $\int_{x}p_{0}c_{i,j}$

Integrating Eq. (155) with respect to $x$, inserting Eq. (156) and noting $p(y)\equiv G(y)$ and $1\equiv H_{0}(y)$ leads to

[TABLE]

Since the functions $H_{i}$ are independent, it first follows

[TABLE]

where $c_{0,0}\equiv 1$, implying $\int_{x}p_{0}c_{0,0}=1$, has been used in the first equation. As these equations must hold for arbitrary $\theta$, it further follows

[TABLE]

Appendix G Autocorrelation of $X(t)$

With the autocorrelation of $X(t)$ given by

[TABLE]

one first finds

[TABLE]

Using $p(x^{\prime},\tau;x,0)\!=\!p(x,0)p(x^{\prime},\tau|x,0)$ then yields the connection to $M^{(1)}$,

[TABLE]

With Eq. (80) one thus finds [up to order $O(\varepsilon^{3})$]

[TABLE]

with unknown but constant coefficients $\lambda_{i}$.