finite-range memory process with zero covariance function after the first
$N+R-1$ values. Using Proposition 4.1.2,
$\boldsymbol{\eta}_{k}\cdot\boldsymbol{x}$ is an MA($N+R-1$) process. This
concludes the proof. ∎

Using Eq. 4.30 the AR($N$) coefficients can be given explicitly

$\phi_{k}=(-1)^{k+1}\sum_{D_{k}}\mathrm{e}^{\Delta t\sum_{i\in
D_{k}}\nu_{i}},$ (4.35)

where $\nu_{i}$ are eigenvalues of the evolution matrix $A$ and $D_{k}$
denotes the family of all $k$-element subsets of the set $\\{1,2,\ldots,N\\}$.
These numbers do not depend on the direction of the projection
$\boldsymbol{x}$ or the exact form of $(\boldsymbol{\eta}_{k})_{k}$’s memory.
The global, deterministic dynamics of the whole state $\boldsymbol{X}$ is
reflected in the AR coefficients. The random dynamics and specific choice of
the observed quantity affects only the MA coefficients. These can be
calculated in two steps, first calculating the covariance of the noise from
Eq. 4.33 and then solving the formula for the covariance 4.16. The solutions
are not unique, see the remark after Proposition 4.3.2. For systems with a
small number of degrees of freedom these solutions are explicit, but often
complicated.

The above theorem can be directly used to study the solutions of higher-order
stochastic equations

######  Proposition 4.2.1.

Let $X$ be a stationary solution of the linear stochastic equation

$a_{q}\frac{\mathrm{d}^{q}}{\mathrm{d}t^{q}}X+a_{q-1}\frac{\mathrm{d}^{q-1}}{\mathrm{d}t^{q-1}}X+\ldots
a_{0}X=f,\quad a_{q}\neq 0$ (4.36)

where $f$ is a Gaussian stochastic stationary force. Then:

  * i)

If $f$ has a covariance with range $R\Delta t$ then the sample $(X_{k})_{k}$
is ARMA($q,q+R-1$).

  * ii)

If $f$ has a covariance which is a sum of $K$ exponential functions then the
sample $(X_{k})_{k}$ is ARMA($q+K,q+K-1$).

###### Proof.

The above differential equation can be expressed as a $q$ dimensional system
of equations with $X_{1}=X$, $X_{2}=\mathrm{d}/\mathrm{d}tX,\ldots$

$\displaystyle\mathrm{d}X_{q-1}$
$\displaystyle=\mathrm{d}t\left(\frac{a_{q-1}}{q_{q}}X_{q-2}+\ldots\frac{a_{0}}{a_{q}}X_{1}\right)+\frac{1}{a_{q}}\mathrm{d}tf$
$\displaystyle\mathrm{d}X_{q-2}$ $\displaystyle=\mathrm{d}tX_{q-1}$
$\displaystyle\vdots$ $\displaystyle\mathrm{d}X_{2}$
$\displaystyle=\mathrm{d}tX_{1},$ (4.37)

Point i) follows immediately. As for ii), a process $f$ with a covariance
function given by $K$ exponents can be decomposed as a sum of $K$ Ornstein-
Uhlenbeck processes $Y_{1},\ldots,Y_{K}$. Therefore the whole systems can be
described by $q+K$ equations

$\displaystyle\mathrm{d}X_{q-1}$
$\displaystyle=\mathrm{d}t\left(\frac{a_{q-1}}{q_{q}}X_{q-2}+\ldots\frac{a_{0}}{a_{q}}X_{1}+\frac{1}{a_{q}}Y_{1}+\ldots+\frac{1}{a_{q}}Y_{K}\right)$
$\displaystyle\mathrm{d}X_{q-2}$ $\displaystyle=\mathrm{d}tX_{q-1}$
$\displaystyle\vdots$ $\displaystyle\mathrm{d}X_{2}$
$\displaystyle=\mathrm{d}tX_{1}$ $\displaystyle\mathrm{d}Y_{1}$
$\displaystyle=-\mathrm{d}t\lambda_{1}Y_{1}+D_{1}\mathrm{d}B_{1}$
$\displaystyle\vdots$ $\displaystyle\mathrm{d}Y_{K}$
$\displaystyle=-\mathrm{d}t\lambda_{K}Y_{K}+D_{K}\mathrm{d}B_{K}.$ (4.38)

Because all noises $\mathrm{d}B_{k}$ are white, the result follows. ∎

###  4.3 Behaviour of sampled Langevin equations

The results from the last section have immediate applications. Many physical
systems are commonly modelled using stochastic linear equations, the most
important ones related to the Langevin and generalized Langevin equation. Few
important examples are given below.

The classical Langevin equation. The case when the driven noise is white was
already considered in Proposition 4.1.1, but different memory models can also
be studied. One of the most important ones is when the driving force is the
Ornstein-Uhlenbeck process, i.e. the coordinate fulfils a set of equations

$\displaystyle\mathrm{d}X$
$\displaystyle=-\mathrm{d}t\lambda_{1}X+\mathrm{d}tF,$
$\displaystyle\mathrm{d}F$
$\displaystyle=-\mathrm{d}t\lambda_{2}F+\mathrm{d}B.$ (4.39)

The eigenvalues of the evolution matrix are very simple in this case

$\nu_{1}=-\lambda_{1},\quad\nu_{2}=-\lambda_{2}.$ (4.40)

The sampled system is ARMA(2,1) with AR coefficients equal to

$\displaystyle\phi_{1}$ $\displaystyle=\mathrm{e}^{-\Delta
t\lambda_{1}}+\mathrm{e}^{-\Delta t\lambda_{2}},$ $\displaystyle\phi_{2}$
$\displaystyle=-\mathrm{e}^{-\Delta t(\lambda_{2}+\lambda_{1})}.$ (4.41)

The inverse relations, useful for estimation, are

$\lambda_{1,2}=-\frac{1}{\Delta
t}\ln\left(\frac{1}{2}\left(\phi_{1}\pm\sqrt{\phi_{1}^{2}+4\phi_{2}}\right)\right).$
(4.42)

An interesting property, visible from the above formulas, is that
$\lambda_{1}$ and $\lambda_{2}$ are not physically distinguishable from AR
coefficients. It is a consequence of the fact that they depend only on the
content of the spectrum, not the order of eigenvalues. To break this symmetry
the additional information contained in MA(1) coefficient must be used.

The full form of classical Langevin equation. In the short time scales the
motion of particle diffusing in water can be modelled by equation

$m\frac{\mathrm{d}^{2}}{\mathrm{d}t^{2}}X=-\kappa
X-\beta\frac{\mathrm{d}}{\mathrm{d}t}X+F,$ (4.43)

governed by the finite-range memory force $F$. It was mentioned as a
motivation at the beginning of the last section. The equation for the charge
$Q$ in a linear RLC circuit with stochastic electromotive force $\mathcal{E}$
has the same form

$L\frac{\mathrm{d}^{2}}{\mathrm{d}t^{2}}Q=-\frac{1}{C}Q-R\frac{\mathrm{d}}{\mathrm{d}t}Q+\mathcal{E}.$
(4.44)

The eigenvalues of the corresponding evolution matrix for both equations, up
to a change of letters, are

$\nu_{1,2}=-\frac{\beta}{2m}\pm\sqrt{\left(\frac{\beta}{2m}\right)^{2}-\frac{\kappa}{m}}.$
(4.45)

The sampled system is ARMA(2,1) with AR coefficients equal to

$\displaystyle\phi_{1}$ $\displaystyle=2\exp\left(-\Delta
t\frac{\beta}{2m}\right)\cosh\left(\Delta
t\sqrt{\left(\frac{\beta}{2m}\right)^{2}-\frac{\kappa}{m}}\right),$
$\displaystyle\phi_{2}$ $\displaystyle=-\exp\left(-\Delta
t\frac{\beta}{m}\right).$ (4.46)

The inverse relations are

$\displaystyle\frac{\beta}{m}$ $\displaystyle=-\frac{1}{\Delta
t}\ln(-\phi_{2}),$ $\displaystyle\frac{\kappa}{m}$
$\displaystyle=-\frac{1}{2\Delta t}\ln(-\phi_{2})-\frac{1}{\Delta
t^{2}}\left(\cosh^{-1}\left(\frac{\phi_{1}}{2\phi_{2}^{2}}\right)\right)^{2}.$
(4.47)

The mass $m$ or inductance $L$ can be determined exclusively using the
variance of the process, as the MA(1) coefficient also depends only on
$\beta/m$ and $\kappa/m$. It is often hard to determine from the experimental
data, because it depends on the units and scale of the system.

RLC circuit with leakage.  The loss of charge is determined by the constant
$G$. In more realistic model the charge may escape from the system. In this
situation the equation for the current $I$ changes and the whole system is
described by

$\displaystyle\mathrm{d}Q$
$\displaystyle=\mathrm{d}tI-\mathrm{d}t\frac{G}{C}Q$
$\displaystyle\mathrm{d}I$
$\displaystyle=-\mathrm{d}t\frac{1}{LC}Q-\mathrm{d}t\frac{R}{L}I+\mathrm{d}\mathcal{E}\frac{1}{L}$
(4.48)

The eigenvalues have a similar form as for the simple RLC circuit

$\nu_{1,2}=-\left(\frac{G}{C}+\frac{R}{L}\right)\pm\sqrt{\left(\frac{1}{2}\left(\frac{G}{C}+\frac{R}{L}\right)\right)^{2}-\frac{1}{CL}},$
(4.49)

therefore the formulas 4.3 and 4.3 for the AR(2) coefficients can be used
after a simple substitution $\frac{R}{L}+\frac{G}{C}\to\frac{\beta}{m}$ and
$\frac{1}{LC}\to\frac{\kappa}{m}$. These two numbers are not sufficient to
determine $R,C$ and $G$, the additional information contained in the MA(1)
process must be used during the estimation.

Generalised Langevin equation with exponential-type kernel. Precisely, the
kernel and the covariance function of the force can be a finite sum of
exponents multiplied by polynomials and sine or cosine functions. This example
is a realisation of our motivation stated at the beginning of Section 4.1.
However, the proper procedure that translates the sampled solution of the GLE
into an ARMA time series is much more complicated than the simple Euler
discretisation would suggest.

To see this we need to show that the solution of the GLE as stated above can
be described in terms of a system of linear equations with constants
coefficients, which will allow us to use Theorem 4.2.2.

######  Theorem 4.3.1.

Let $X$ be a stationary solution of the GLE

$\frac{\mathrm{d}}{\mathrm{d}t}X=-\int_{-\infty}^{t}\\!\\!\\!\mathrm{d}\tau\
X(\tau)K(t-\tau)+F,$ (4.50)

where the kernel $K$ and the covariance function $r_{F}$ are finite sums of
exponentials, which can be multiplied by polynomials and sine or cosine
functions. In such case $X$ is equivalent to
$\boldsymbol{Y}\cdot\boldsymbol{q}$, where $X$ is a solution of a system of
linear stochastic equations with constant coefficients, and $\boldsymbol{q}$
come fixed vector.

Proof. The Fourier transform of $\sin(vt)\mathrm{e}^{-\lambda t}$ on
$\mathbb{R}_{+}$ is $v\big{(}v^{2}+(\lambda-i\omega)^{2}\big{)}^{-1}$. The
Fourier transforms of $K$ and $r_{F}$ are sums of derivatives and polynomial
times this quantity, so they are rational functions. From our considerations
in Section 2.2 and 3.4 we know that the power spectral density of the
stationary solution is

$s_{X}(\omega)=\frac{s_{F}}{|\mathrm{i}\omega+\widehat{K}|^{2}}=\frac{P(\omega)}{Q(\omega)}$
(4.51)

i.e., as a ratio of two rational functions, it is also a rational function and
it may be assumed that it is a ratio of the polynomials $P$ and $Q$ with no
common roots. Moreover, this function must be even, like any power spectral
density. It is even if and only if $P$ and $Q$ are even, because
$P(\omega)Q(-\omega)=P(-\omega)Q(\omega)$ implies that $P(\omega)$ and
$P(-\omega)$ have the same roots.

It may not be immediately clear, but both $P$ and $Q$ can be taken to be
square modulus of some other complex valued polynomials. This property holds
because they are positive and even, so they have only complex roots that come
in conjugate pairs. They are products of terms like $|(\omega-z_{k})(\omega-
z_{k}^{*})|$, and as $\omega$ is real, any such term is equal to $|\omega-
z_{k}|^{2}=|\omega-z_{k}^{*}|^{2}$. Moreover, in this decomposition we can
only choose roots $z_{k}$ with positive imaginary part.

Now, let us consider $P$. We can write
$P(\omega)=\big{|}\sum_{k=0}^{n}(a_{k}+\mathrm{i}b_{k})\omega^{k}\big{|}^{2},a_{k},b_{k}\in\mathbb{R}$.
Therefore

$P(\omega)=\big{(}a_{n}\omega^{n}+a_{n-1}\omega^{n-1}+\ldots+a_{0}\big{)}^{2}+\big{(}b_{n}\omega^{n}+b_{n-1}\omega^{n-1}+\ldots+b_{0}\big{)}^{2}$
(4.52)

and because $P$ is even, the above formula is a function of $\omega^{2}$. The
left and right parentheses also must be a polynomials of $\omega^{2}$ which is
only possible when any of them contains only odd or even powers of $\omega$
inside. In other words $a_{2k}=0$ and $b_{2k+1}=0$ or $a_{2k+1}=0$ and
$b_{2k}=0$. So $P$ has the form

$P(\omega)=\left|\sum_{k=0}^{n}p_{k}(\mathrm{i}\omega)^{k}\right|^{2},\quad\tilde{p}_{2k}=(-1)^{k}a_{2k},p_{2k+1}=(-1)^{k}b_{2k+1},p_{k}\in\mathbb{R}.$
(4.53)

The same goes for
$Q(\omega)=\big{|}\sum_{j=0}^{m}q_{j}(\mathrm{i}\omega)^{j}\big{|}^{2},q_{j}\in\mathbb{R}$.
Consider a process $X$ defined by two equations

$\displaystyle
p_{n}\frac{\mathrm{d}^{n}}{\mathrm{d}t^{n}}Y+p_{n-1}\frac{\mathrm{d}^{n-1}}{\mathrm{d}t^{n-1}}Y+\ldots+p_{0}Y=\frac{\mathrm{d}}{\mathrm{d}t}B,$
(4.54) $\displaystyle
X=q_{m}\frac{\mathrm{d}^{m}}{\mathrm{d}t^{m}}Y+q_{m-1}\frac{\mathrm{d}^{n-1}}{\mathrm{d}t^{n-1}}Y+\ldots+q_{0}Y.$
(4.55)

We have chosen the polynomial
$\omega\mapsto\sum_{k=0}^{n}p_{k}(\mathrm{i}\omega)^{k}$ such that it has only
roots with positive imaginary parts. This implies that the polynomial
$s\mapsto\sum_{k=0}^{n}p_{k}s^{k}$ has roots with negative real parts. It
guarantees that Eq. 4.54 has a stationary solution. The process $X$ is well-
defined as $m<n$ and also stationary. It has exactly the psd 4.51 and has the
same distribution as the solution of given GLE.

Equations 4.54 and 4.55 are equivalent to a vector equation for
$\boldsymbol{Y}$ made of $Y$ and its $n-2$ derivatives, additionally
$X=\boldsymbol{Y}\cdot\boldsymbol{q}$ where
$\boldsymbol{q}=[0,0,\ldots,0,q_{m},q_{m-1},\ldots,q_{0}]$.

$\blacksquare$

The hard part of the proposition above is proving that the reduction of the
GLE to a system of linear equations is always possible, given the assumptions.
The procedure itself is not so complex, it requires mainly calculating
polynomials $P,Q$ and their roots.

As an example, let us consider the GLE with exponential kernel:
$K(t)=b^{2}\mathrm{e}^{-2at}$ for which the fluctuation-dissipation theorem
holds. We have already provided some basic information about this model
following Eq. 3.73. Now we consider its Fourier space representation,

$\widehat{K}(\omega)=\frac{b^{2}}{2a-\mathrm{i}\omega},\quad
s_{F}(\omega)=\widehat{r}_{F}(\omega)=\frac{4ab^{2}}{4a^{2}+\omega^{2}},$
(4.56)

consequently

$s_{X}(\omega)=\frac{4ab^{2}}{4a^{2}+\omega^{2}}\frac{|2a-\mathrm{i}\omega|^{2}}{|\mathrm{i}\omega(2a-\mathrm{i}\omega)+b^{2}|^{2}}=\frac{4ab^{2}}{|\omega^{2}+\mathrm{i}2a\omega+b^{2}|^{2}}.$
(4.57)

The roots of polynomial $\omega^{2}+2\mathrm{i}a\omega+b^{2}$ do not both have
positive imaginary part

$\omega^{2}+2\mathrm{i}a\omega+b^{2}=\left(\omega-\mathrm{i}\big{(}\lambda+\sqrt{\lambda^{2}+b^{2}}\big{)}\right)\left(\omega-\mathrm{i}\big{(}\lambda-\sqrt{\lambda^{2}+b^{2}}\big{)}\right).$
(4.58)

The first root has positive imaginary part, the second one negative, so it
should be conjugated,

$\left(\omega-\mathrm{i}\big{(}\lambda+\sqrt{\lambda^{2}+b^{2}}\big{)}\right)\left(\omega-\mathrm{i}\big{(}-\lambda+\sqrt{\lambda^{2}+b^{2}}\big{)}\right)=\omega^{2}-2\mathrm{i}\sqrt{\lambda^{2}+b^{2}}\omega-b^{2}.$
(4.59)

The solution of the GLE is equivalent to a solution of the stochastic
differential equation

$\frac{\mathrm{d}^{2}}{\mathrm{d}t^{2}}X+2\sqrt{\lambda^{2}+b^{2}}\frac{\mathrm{d}}{\mathrm{d}t}X+b^{2}X=2b\sqrt{a}\frac{\mathrm{d}}{\mathrm{d}t}B$
(4.60)

which is a full form of the classical Langevin equation with stiffness
$\kappa=b^{2}$ and friction coefficient $\beta=2\sqrt{a^{2}+b^{2}}$. The
sampled solution of this GLE is an ARMA(2,1) process with AR coefficients

$\displaystyle\phi_{1}$ $\displaystyle=2\exp\left(-\Delta
t\sqrt{a^{2}+b^{2}}\right)\cosh\left(\Delta ta\right),$
$\displaystyle\phi_{2}$ $\displaystyle=-\exp\left(-2\Delta
t\sqrt{a^{2}+b^{2}}\right).$ (4.61)

As we see during a statistical analysis $a$ can be estimated using
$\phi_{1}^{2}/\phi_{2}$.

Studying models, such as the examples given above, can be performed in both
continuous- and discrete-time. To some degree both approaches are equivalent.
For example, the covariance function of the sampled series $(X_{k})_{k}$ is
simply the covariance of the continuous-time $(X(t))_{t}$ at points $\Delta
tk$, i.e. $(r_{X}(\Delta tk))_{k}$. But there are also differences.

One method, which has no variant in the continuous-time case, is studying a
partial autocorrelation function (pacf) $\rho_{X}$. It is the memory function
measuring a dependence between $X_{i}$ and $X_{i+k}$ with influences from in-
between $X_{j},i<j<i+k$ removed. This removal is performed by subtraction of
the projection onto the subspace of the variables ${X_{i+1},\ldots,X_{i+k}}$.
Explicitly:

######  Definition 4.3.1.

The partial autocorrelation function of the stationary time series $X$,
$\rho_{X}(k)$ is given by formula

$\rho_{X}(k)\mathrel{\mathop{:}}=\operatorname*{{corr}}(X_{i}-\mathbb{E}[X_{i}|\\{X_{j}\\}_{j=i+1}^{i+k-1}],X_{i+k}-\mathbb{E}[X_{i+k}|\\{X_{j}\\}_{j=i+1}^{i+k-1}]),$
(4.62)

where we additionally use the convention that
$\mathbb{E}[\boldsymbol{\cdot}|\varnothing]=0$.

The partial autocorrelation is a useful quantity, because by removal of these
projections we get rid of the influence of the time evolution between instants
$i$ and $j$. This procedure helps to clarify the analysis of the memory. The
pacf is a measure of the direct dependence between the values of the process
and is well-suited to study AR time series.

######  Proposition 4.3.1.

The partial autocorrelation of an AR($p$) time series is zero for lag equal or
greater than $p$.

###### Proof.

The representation

$X_{k}=\phi_{1}X_{k-1}+\ldots+\phi_{p}X_{k-p}+\xi_{k}$ (4.63)

decomposes $X_{k}$ as a sum of $\xi_{k}$, which is independent from the
$\mathrm{span}\\{X_{k-1},X_{k-2},\ldots\\}$ and
$\phi_{1}X_{k-1}+\ldots+\phi_{p}X_{k-p}$, which is clearly a variable in
$\mathrm{span}\\{X_{k-1},X_{k-2},\ldots,X_{k-p}\\}$. So, for $r\geq p$

$\mathbb{E}[X_{k}|X_{k-1},\ldots,X_{k-r}]=\phi_{1}X_{k-1}+\ldots+\phi_{p}X_{k-p}$
(4.64)

and

$\rho_{X}(r)=\operatorname*{{corr}}(\xi_{k},X_{k-r}-\mathbb{E}[X_{k-r}|X_{k-1},\ldots,X_{k-r+1}])=0$
(4.65)

due to the aforementioned independence. ∎

For the AR(1) process the one significant non-zero value of partial
autocorrelation is easy to calculate,
$\rho_{X}(1)=\operatorname*{{corr}}(\phi_{1}X_{k-1}+\xi_{k},X_{k-1})=\phi_{1}$,
so

$\rho_{X}(k)=\begin{cases}1,&k=0;\\\ \phi_{1},&k=1;\\\ 0,&k>1.\end{cases}$
(4.66)

In the general case the partial autocorrelation can be estimated using the
Yule-Walker equations, which relate it to the covariance function (Brockwell
and Davis, 2006). This set of equations can be easily solved, at least
numerically.

Apart from the partial autocorrelation function, there is a significant
difference in the spectral form of the continuous and discrete-time processes.
The latter, instead of the whole $\mathbb{R}$ has a spectral measure defined
on a torus (a natural dual space), or, equivalently, on the interval
$[-\pi,\pi]$, see Theorem 2.2.2. Indeed, a short calculation yields

$\displaystyle X_{k}$ $\displaystyle=X(\Delta
tk)=\int_{\mathbb{R}}\sigma(\mathrm{d}\omega)\
\mathrm{e}^{\mathrm{i}\omega\Delta t}=\int_{-\frac{\pi}{\Delta
t}}^{\frac{\pi}{\Delta
t}}\left(\sum_{j=-\infty}^{\infty}\sigma\left(\mathrm{d}\omega+\frac{2\pi}{\Delta
t}k\right)\right)\ \mathrm{e}^{\mathrm{i}\omega\Delta tk}$
$\displaystyle=\frac{1}{\Delta
t}\int_{-\pi}^{\pi}\left(\sum_{j=-\infty}^{\infty}\sigma\left(\frac{1}{\Delta
t}(\mathrm{d}\omega+2\pi k)\right)\right)\ \mathrm{e}^{\mathrm{i}\omega k},$
(4.67)

which relates the spectral measure of a process and the corresponding sampled
time series. Only the discrete variant can be estimated from the data, but, as
one can see, the full information about the process can be regained in the
infill limit $\Delta t\to 0$. Equation 4.3 can be used to calculate the
spectral measure of the sampled process, but in most of the cases the infinite
sum within is hard to study analytically. For an ARMA processes there exists a
much stronger result.

######  Proposition 4.3.2.

The power spectral density $s$ of an ARMA($p,q$) process is

$s(\omega)=c\frac{\left|1+\sum_{k=1}^{q}\theta_{k}\mathrm{e}^{-\mathrm{i}\omega
kt}\right|^{2}}{\left|1-\sum_{k=1}^{p}\phi_{k}\mathrm{e}^{-\mathrm{i}\omega
kt}\right|^{2}},$ (4.68)

where $c$ is a constant spectral density of $(\xi_{k})_{k}$.

###### Proof.

This proposition is similar to results for linear filters in Section 2.2.
Namely, we use the harmonic representation for the series $(X_{k})_{k}$, as
stated in Theorem 2.2.2

$X_{k}=\int_{-\pi}^{\pi}\mathrm{d}S_{X}(\omega)\ \mathrm{e}^{\mathrm{i}\omega
k}.$ (4.69)

From this formula it is clear that the process $(X_{k-1})_{k}$ has the
spectral process $\mathrm{d}S_{X}(\omega)\mathrm{e}^{-\mathrm{i}\omega}$. A
similar shift occurs for $(\xi_{k})_{k}$, so the following relation holds

$\int_{-\pi}^{\pi}\mathrm{d}S_{X}(\omega)\left(1-\sum_{k=1}^{p}\phi_{k}\mathrm{e}^{-\mathrm{i}\omega
kt}\right)\mathrm{e}^{\mathrm{i}\omega
t}=\int_{-\pi}^{\pi}\mathrm{d}S_{\xi}(\omega)\left(1+\sum_{k=1}^{q}\theta_{k}\mathrm{e}^{\mathrm{i}\omega
kt}\right)\mathrm{e}^{\mathrm{i}\omega t}.$ (4.70)

The uniqueness of the harmonic representation guarantees that

$\mathrm{d}S_{X}(\omega)=\mathrm{d}S_{\xi}(\omega)\frac{1+\sum_{k=1}^{q}\theta_{k}\mathrm{e}^{-\mathrm{i}\omega
kt}}{1-\sum_{k=1}^{p}\phi_{k}\mathrm{e}^{-\mathrm{i}\omega kt}}.$ (4.71)

The power spectral density is the variance of $\mathrm{d}S_{X}$, so the result
holds. ∎

A consequence of this proposition is that the coefficients of the MA processes
are not uniquely determined, which complicates their physical interpretation.
A distribution of the process is uniquely determined by its spectral measure,
which for MA process is
$\left|1+\sum_{k=1}^{q}\theta_{k}\mathrm{e}^{-\mathrm{i}\omega
kt}\right|^{2}$. It is a polynomial of the argument
$\mathrm{e}^{-\mathrm{i}\omega}$, which has $q$, not necessarily distinct,
zeros $z_{k}$. It can be factorised into a product of $q$ terms
$\left|z_{k}-\mathrm{e}^{-\mathrm{i}\omega}\right|^{2}$. However, because of
the modulus, if the zero $z_{k}$ is real, the corresponding factor can be
transformed as

$\left|z_{k}-\mathrm{e}^{-\mathrm{i}\omega}\right|=\left|1-z_{k}\mathrm{e}^{\mathrm{i}\omega}\right|=|z_{k}|\left|z_{k}^{-1}-\mathrm{e}^{-\mathrm{i}\omega}\right|.$
(4.72)

Note that factor $|z_{k}|$ only changes the variance of the resulting process.
For complex zeros this procedure would result in a complex polynomial and a
complex MA process if $z_{k}$ and its complex conjugate $z_{k}^{*}$ are not
inversed at the same time. Because we consider only real-valued MA processes,
the number of possibilities is $2^{q^{\prime}}$, where $q^{\prime}$ is a
number of $z_{k}$ which are real plus half of the number of complex ones.

For a real-valued MA(1) process there are always exactly two options:
$\theta_{1}$ and $\theta_{1}^{-1}$. The corresponding processes differ only by
the variance, which changes by the ratio $\theta_{1}^{2}(1+\theta_{1}^{2})$.

###  4.4 Application to modelling of optical tweezers

An example of a system in which our approach is especially viable are optical
tweezers. These are a versatile tool allowing to manipulate micrometer-sized
particles in liquids (Neuman and Block, 2004) non-invasively and to measure
forces even on the picoNewton scale. For this reason they have a wide range of
applications in many fields of biology and soft condensed matter physics
including, e.g. stretching of DNA (Wang et al., 1997) and other polymers
(Smith et al., 1996), molecular motors (Arsenault et al., 2009) research or
analysis of colloidal suspensions (Meyer et al., 2006).

The most classical form of the measurements related to optical tweezers are
series of of positions of a bead trapped in the tweezers. Observations
performed using CCD camera are especially useful, because the camera provides
a lot of diverse information and allows for the tracking of many objects
simultaneously. The data that we use was measured using this setup by Sławomir
Drobczyński at Wrocław University of Science and Technology. As we will show
in our analysis below, in the case of high-frequency measurements the camera
has a significant influence on the measurements, which is possible to analyse
using discrete-time methods.

According to the classical Einstein theory of diffusion, a trajectory
$(X(t))_{t\in\mathbb{R}}$ of a colloidal particle trapped in a viscous fluid
by optical tweezers is a solution of the equation of the force balance

$0={F}_{S}+{F}_{O}+{F}_{T},$ (4.73)

where

  * -

$F_{S}=-\beta\frac{\mathrm{d}X}{\mathrm{d}t}$ is the Stokes force (friction of
liquid) acting on the spherical bead with radius $r$ within the liquid with
viscosity $\eta$, $\beta=6\pi\eta r$;

  * -

$F_{O}=-\kappa X$ is the force caused by the optical tweezers, we use an
harmonic approximation, i.e., assume that the potential is harmonic with
stiffness $\kappa$;

  * -

$F_{T}=\sqrt{k_{B}\mathcal{T}\beta}{\mathrm{d}B}/{\mathrm{d}t}$ is the thermal
force: it models the exchange of momenta with particles of the liquid.

Substitution of the explicit formulas assures that the above force balance
condition is equivalent to the Langevin equation

$\mathrm{d}X=-\mathrm{d}t\lambda X+D\mathrm{d}B,$ (4.74)

where $\lambda=\kappa/\beta$ and $D=\sqrt{k_{B}\mathcal{T}/\beta}.$ So, the
stationary solution is the Ornstein-Uhlenbeck process, and the sampled time
series is an AR(1) process. The AR(1) coefficient $\phi_{1}$ can be estimated
using the autoregressive nature of the process. One can use the classical
regression with the sequence of $(X_{k})_{k}$ taken as a variable y and
shifted sequence of $(X_{k-1})_{k}$ taken as a variable x. The well-known
least-squares estimator of the slope in this case reads

${\phi}_{1}^{\text{est}}=\frac{\sum_{k}X_{k}X_{k-1}}{\sum_{k}X_{k}^{2}}.$
(4.75)

However, using a partial autocorrelation is a more robust method, because the
estimated $\rho_{X}(1)$ is an estimate of $\phi_{1}$ (Eq. 4.66) and checking
if $\rho_{X}(k)\approx 0$ for $k>1$ is a form of the test of validity of the
model.

However, the pure physical model of the optical tweezers is not sufficient to
analyse typical recordings because it does not take into account the
distortions which are caused by the video camera. The most basic and well-
researched type of distortion is an addition of the noise generated by the
camera’s CMOS matrix. In this case we do not observe $X_{k}$, but the series
$X^{\prime}_{k}=X_{k}+w_{k}$, where $w_{k}$ is a noise independent from the
values of $X_{k}$. The independence of $X_{k}$ and $w_{k}$ implies that the
resulting covariance is the sum $r_{X^{\prime}}=r_{X}+r_{w}$, the same for the
spectral measure $\sigma_{X^{\prime}}=\sigma_{X}+\sigma_{w}$. In order to
understand the structure of $(X^{\prime}_{k})_{l}$ let us look what linear
filter governs its evolution

$\displaystyle X^{\prime}_{k}$
$\displaystyle=X_{k}+w_{k}=\phi_{1}X_{k-1}+\xi_{k}+w_{k}$
$\displaystyle=\phi_{1}X^{\prime}_{k-1}+\xi_{k}+w_{k}-\phi_{1}w_{k-1}=\phi_{1}X^{\prime}_{k-1}+\xi_{k}^{\prime}.$
(4.76)

Now, we take the common assumption that the additive noise $(w_{k}))k$ is
Gaussian and white, with a variance $s^{2}$. Then the new effective noise
$\xi^{\prime}_{k}=\xi_{k}+w_{k}-\phi_{1}w_{k-1}$ has the covariance function

$r_{\xi^{\prime}}(h)=\begin{cases}\frac{D^{2}}{2\lambda}\left(\phi_{1}^{-2}-1\right)+\left(1+\phi_{1}^{2}\right)s^{2},&h=0;\\\
-\phi_{1}s^{2},&h=1;\\\ 0,&h>1\end{cases}$ (4.77)

As we see it is a time series with a short-term negative memory. As stated
before, any such process can be expressed as a linear filter acting on a
Gaussian white noise, in this case it is an MA(1) time series

$\xi^{\prime}_{k}=\eta_{k}-a\eta_{k-1},$ (4.78)

where $(\eta_{k})_{k}$ is a white noise series with variance $v^{2}$. To
present the constant $a$ in a simple form it is convenient to introduce a
ratio

$\mathfrak{r}=-\frac{r_{\xi^{\prime}}(0)}{r_{\xi^{\prime}}(1)},$ (4.79)

which is always greater or equal than two, because $\phi_{1}>0$ implies
$\mathfrak{r}\geq\phi_{1}+\phi_{1}^{-1}\geq 2$. Solving the quadratic equation
obtained comparing the covariances, we get only one solution with $|a|<1$,

$a=\frac{\mathfrak{r}-\sqrt{\mathfrak{r}^{2}-4}}{2},\ \ \
v^{2}=\frac{2r_{\xi^{\prime}}(0)}{\mathfrak{r}(\mathfrak{r}-\sqrt{\mathfrak{r}^{2}-4})}.$
(4.80)

We may consider the series $\eta_{k}$ as the noises $\xi_{k}$ and $w_{k}$
mixed and orthogonalised. Through this calculation we have shown that the
observed $X^{\prime}_{k}$ fulfils

$X^{\prime}_{k}-\phi_{1}X^{\prime}_{k-1}=\eta_{k}-a\eta_{k-1}.$ (4.81)

This is a specific ARMA(1,1) time series.

Let us once again turn back to the influence of the camera. Additive white
noise is not the only common type of distortion. The cameras have a tendency
to blur the images. Especially when the frequency of sampling is high and the
exposure time is long, the CMOS matrix may not have sufficient time to fully
refresh. As a result, the remainder of a previous frame is left in the image
of the next one. Instead of $X_{k}$ we observe $\tilde{X}_{k}=X_{k}+bX_{k-1}$,
$0<b<1$. It is also possible that $\tilde{X}_{k}$ depends on even older
images, but the intensity of this effect should decay exponentially. Our
measurements show that the observed $b$ is not very large, so we will neglect
this effect.

So, let us denote a short-time linear filter $M=1+b\mathcal{S}_{-1}$. The
operators $M$ and $L=1-\phi_{1}\mathcal{S}_{-1}$ commute and

$L\tilde{X}_{k}=LMX_{k}=MLX_{k}=M\xi_{k}$ (4.82)

or, equivalently,

$\tilde{X}_{k}-\phi_{1}\tilde{X}_{k-1}=\xi_{k}+b\xi_{k-1}.$ (4.83)

Once again the observed process is ARMA(1,1), but this time with the positive
MA(1) coefficient equal to $b$. We have come to the conclusion that both
possible distortions caused by the digital camera cause the observed series to
fulfil the ARMA(1,1) equation, the difference being the sign of the MA(1)
part.

Let us denote the observed series by $\mathcal{X}_{k}$, which may be
$X^{\prime}_{k}$ or $\tilde{X}_{k}$. The direct dependence between
$\mathcal{X}_{k}$ and the noise $\eta_{k}$ can be studied by applying the
inverse operator
$(1-\phi_{k}\mathcal{S}_{-1})^{-1}=\sum_{h=0}^{\infty}\phi_{1}^{h}\mathcal{S}_{-h}$
to both sides of the above equation, yielding

$\mathcal{X}_{k}=\eta_{k}+\left(\phi_{1}+\theta_{1}\right)\sum_{j=1}^{\infty}\phi_{1}^{j-1}\eta_{k-j}.$
(4.84)

The covariance of this series is straightforward to calculate and equal to

$r_{\mathcal{X}}(h)=v^{2}\cdot\begin{cases}1+\frac{(\phi_{1}+\theta_{1})^{2}}{1-\phi_{1}^{2}},&h=0;\\\
\phi_{1}^{|h-1|}(\phi_{1}+\theta_{1})\frac{1+\phi_{1}\theta_{1}}{1-\phi_{1}^{2}},&h\neq
0;\end{cases}$ (4.85)

Hence, the ARMA(1,1) process also exhibits an exponential decay of the
covariance function with one special value at $h=0$. It is a natural insight
for the case $\mathcal{X}_{k}=X^{\prime}_{k}$ when we can simply add the
covariance function of the original process and white noise, but is also true
in a more general situation. It can be interpreted as an effect of distortion
which affects only the shortest one-time-step memory.

Similarly to the previously considered problems, the power spectral density of
ARMA(1,1) process can be obtained using Proposition 4.3.2, and equals

$s_{\mathcal{X}}(\omega)=v^{2}\frac{1+\theta_{1}^{2}+2\theta_{1}\cos(\omega\Delta
t)}{1+\phi_{1}^{2}-2\phi_{1}\cos(\omega\Delta t)}.$ (4.86)

A comparison of the spectra for different processes relevant to the optical
tweezers measurements is shown in Fig. 4.3. One must be careful, because the
psd of the Ornstein-Uhlenbeck process $\propto 1/(\omega^{2}+\lambda^{2})$ is
for low frequencies nearly identical to the truly observed psd disturbed by a
blur described by a small MA coefficient. This is an unfortunate coincidence
which can easily misguide the analysis, creating a false resemblance of the
lack of distortions.

Figure 4.3: Comparison of psd functions for continuous-time Ornstein-Uhlenbeck
process, discrete-time undisturbed and disturbed ARMA(1,1) time series.

To study the approximate behaviour of the partial autocorrelation function we
inverse the operator $1+\theta_{1}\mathcal{S}_{-1}$, expressing
$\mathcal{X}_{k}$ directly as a function of past $\mathcal{X}_{k-j}$,

$\mathcal{X}_{k}=(\phi_{1}+\theta_{1})\sum_{j=1}^{\infty}(-\phi_{1})^{j-1}\mathcal{X}_{k-j}+\eta_{k}.$
(4.87)

It shows that the direct dependence between $\mathcal{X}_{k}$ and
$\mathcal{X}_{k-j}$ is proportional to the constant
$(\phi_{1}+\theta_{1})(-\phi_{1})^{j-1}$; the partial autocorrelation of
$\mathcal{X}$ is a geometrical decay but its rate depends on the MA
coefficient $\theta_{1}$ instead of the AR coefficient $\phi_{1}$ as for the
covariance function. As $\theta_{1}$ is close to one for common sampling
frequencies and $\phi_{1}$ is not large even for highly disturbed data, we
have $\phi_{1}+\theta_{1}>0$ and the sign of $\rho_{\mathcal{X}}(2)$ is the
same as the sign of $-\phi_{1}$. Thus, the estimation of the partial
autocorrelation function immediately determines the presence and type of
distortions in the data. The values of the pacf for larger $h$ have little
practical importance in our case; because of the fast decay these values
quickly become unobservable. The first values of the pacf can be calculated
analytically

$\rho_{\mathcal{X}}(h)=\begin{cases}\frac{(\phi_{1}+\theta_{1})(1+\phi_{1}\theta_{1})}{1+2\phi_{1}\theta_{1}+\phi_{1}^{2}},&h=1;\\\
-\phi_{1}\frac{(\phi_{1}+\theta_{1})(1+\phi_{1}\theta_{1})}{(1+2\phi_{1}\theta_{1}+\phi_{1}^{2})^{2}},&h=2.\end{cases}$
(4.88)

The partial autocorrelation function estimated from the data very closely
follows the described behaviour, an example is shown in Fig. 4.4. From the
rate of the exponential decay we can assess that $\theta_{1}\approx 0.17$.
This is a very good preliminary estimate, even compared to the more precise
methods that we present later.

Figure 4.4: Estimated pacf of the optical tweezers time series; sampling
frequency is 5 [kHz], camera exposure time is 100 [$\mu$s].

Both distortions considered above are very common, so the natural question
appears: what happens if they are present simultaneously? Their joint effect
depends on the sequence in which video camera imperfections distort the
recording. One possibility is that the camera records images disturbed by the
white noise, which situation corresponds to $X^{\prime}_{k}$, and the film is
blurred later, so that we observe $\tilde{X}^{\prime}_{k}=MX^{\prime}_{k}$. In
this case

$\displaystyle L\tilde{X}^{\prime}_{k}$
$\displaystyle=M\xi_{k}^{\prime}=(1+b\mathcal{S}_{-1})(1-a\mathcal{S}_{-1})\eta_{k}$
$\displaystyle=\eta_{k}+(b-a)\eta_{k-1}-ba\eta_{k-2}.$ (4.89)

This is an ARMA(1,2) process. Note that the first MA coefficient
$\theta_{1}=b-a$ has a sign determined by the relative strengths of both
distortions; the second coefficient $\theta_{2}=-ba$ is relatively small and
its presence may be hard to detect.

The situation is different if the camera adds the additive noise to the
already blurred recording, e.g. during processing of the electrical signal
from the CMOS matrix. Then we observe
${\tilde{X}_{k}}^{\prime}=\tilde{X}_{k}+w_{k}$. This process fulfils the
relation

$\displaystyle{\tilde{X}_{k}}^{\prime}$
$\displaystyle=\phi_{1}\tilde{X}_{k}+\xi_{k}+b\xi_{k-1}+w_{k}$
$\displaystyle=\phi_{1}{\tilde{X}}^{\prime}_{k-1}+\xi_{k}+b\xi_{k-1}+w_{k}-\phi_{1}w_{k-1}.$
(4.90)

The new effective noise
$\tilde{\xi}_{k}=\xi_{k}+b\xi_{k-1}+w_{k}-\phi_{1}w_{k-1}$ has the covariance
function

$r_{\tilde{\xi}}(h)=\begin{cases}\left(1+b^{2}\right)\frac{D^{2}}{2\lambda}\left(\phi_{1}^{-2}-1\right)+\left(1+\phi_{1}^{2}\right)s^{2},&h=0;\\\
-\phi_{1}s^{2}+b\frac{D^{2}}{2\lambda}\left(\phi_{1}^{-2}-1\right),&h=1;\\\
0,&h>1.\end{cases}$ (4.91)

This is again an MA(1) process; its $\theta_{1}$ parameter is given by Eqs.
4.79 and 4.80 with $r_{\tilde{\xi}}$ inserted instead of $r_{\xi^{\prime}}$.
The time series $({\tilde{X}_{k}}^{\prime})_{k}$ is ARMA(1,1) but with MA
parameter different from that of $X_{k}^{\prime}$ or $\tilde{X}_{k}$.

The additive noise and blur have opposite influences, which can be seen
looking at the case $h=1$ in 4.91. It is even possible for them to cancel each
other in which case the recording behaves like an undisturbed trajectory of
the stochastic oscillator, i.e. an AR(1) time series. This effect can be
created involuntarily during calibration of the experimental setup. For the
case “blur+noise” (${\tilde{X}}^{\prime}$) there is a possibility for exact
cancelling of the blur and additive noise effects, for the “noise+blur” case
$\tilde{X}^{\prime}$ there will be a small second MA(2) coefficient left equal
to $\theta_{2}=-a^{2}=-b^{2}$. For both cases the only large statistical
change would be an increase of the generating noise variance ($h=0$ in 4.91)
and a subsequent increase of the variance for the process $X$ itself. As
measuring the variance of $X$ is one of the popular methods of estimating the
stiffness of the trap, special care must be taken to avoid the discussed
possibility.

We have come to the conclusion that the presence and type of the experimental
distortions observed in the data can be recognised by the order of the ARMA
time series and the value of the MA coefficient. At the same time, the AR
coefficient $\phi_{1}$ does not depend on the distortions. Hence, the
stiffness coefficient of the tweezers’ harmonic trap can always be accurately
determined through fitting of the ARMA model.

The simplest method to start the analysis is to fit parameters to an estimated
memory function: covariance, partial autocorrelation or power spectral
density. The covariance function was estimated as the sample covariance (e.g.
method acf in R package). The sample pacf was obtained by substituting into
the Yule-Walker equations the values of the sample covariance. This method is
already implemented in the Matlab environment as a parcorr function, and in
the R package as pacf. To estimate the discrete-time psd we used the
periodogram defined as

$I(\omega)=\frac{\Delta
t}{N}\left|\sum_{k=1}^{N}X_{k}\mathrm{e}^{\mathrm{i}\Delta
tk\omega}\right|^{2}.$ (4.92)

The periodogram itself is an ubiased estimator, but is not consistent. For
this reason we take a smoothed periodogram as an estimate, that is the
convolution of the function 4.92 with a smoothing kernel. For a sufficient
choice of smoothing the result has better statistical properties than the pure
periodogram. In Matlab this method is implemented as command periodogram (uses
a rectangular smoothing function) and in the R package as the function
spec.pgram with more options available allowing to regulate the calculation
method and smoothing.

The values of $v^{2},\phi_{1}$ and $\theta_{1}$ obtained from fitting one of
the memory functions are suboptimal estimates. Better results are obtained
using the least-squares method, i.e. looking for the values
${\phi}_{1}^{\text{est}},{\theta}_{1}^{\text{est}}$ for which the sum

$\sum_{k=1}^{N}({\eta}_{k}^{\text{est}})^{2}=\sum_{k=1}^{N}\left(X_{k}-{\phi}_{1}^{\text{est}}X_{k-1}-{\theta}_{1}^{\text{est}}{\eta}_{k-1}^{\text{est}}\right)^{2}$
(4.93)

has a minimal value. Here ${\eta}_{k}^{\text{est}}$ are calculated
recursively; as an initial value ${\eta}_{0}^{\text{est}}$ we take 0. If our
estimates ${\phi}_{1}^{\text{est}}$ and ${\theta}_{1}^{\text{est}}$ are close
to the real values, the estimated series $({\eta}_{k}^{\text{est}})$ is close
to $(\eta_{k})_{k}$ and is similar to white noise. If they are not close, the
variables ${\eta}_{k}^{\text{est}}$ are correlated, which causes the sum of
squares to increase. Minimising the above sum, given initial estimates for
$\phi_{1}$ and $\theta_{1}$ is an accurate and computationally fast method. In
Matlab it is provided by the command estimated and in the R package by arima
or, using maximum likelihood method, in the library FitARMA. The estimate of
the psd obtained from the data and the psd of the fitted ARMA(1,1) process is
shown in Fig. 4.5. This figure proves that the least-squares fit of an ARMA
model provides estimates which fully explain the spectral properties of the
measurements.

Figure 4.5: Smoothed periodogram (blue thin line) of the data and the dpsd of
the ARMA(1,1) process (thick red line), ${\phi}^{\text{est}}_{1}=0.9580$,
${\theta}_{1}^{\text{est}}=0.0472$.

We have experimentally studied the influence of the digital camera for eight
different sampling frequencies and seven different exposure times. The
results, obtained using least-squares fitting, are shown in Fig 4.6. The shape
of this plot is very natural to explain. For short exposure times not much
light goes into camera and additive noise is an important factor compared to
the intensity of the undisturbed image. For longer exposure times there is
more light, but also less time to refresh the CMOS matrix between subsequent
frames, which causes blur to dominate. Increasing sampling frequencies further
shortens time when the CMOS matrix can refresh and strengthens this effect.
The level for which the MA coefficient $\theta_{1}\approx 0$ corresponds to
the situation when these two effects cancel each other, leading to the false
illusion of the distortion absence.

The MA coefficient depends on the sampling frequency $f$ (in kHz) and the
exposure time $t_{e}$ in a regular manner, approximately linearly in the
studied range, so that even a simple linear approximation of the obtained
dependence: $\theta_{1}=0.0075f+0.001t_{e}-0.0798$ provides a good
interpolation. However, our measurements for larger frequencies suggest that
the studied effect becomes non-linear in the larger range.

During these measurements we also checked the values of the estimated AR(1)
coefficient. During few subsequent measurements we obtained estimates
${\phi}_{1}^{\text{est}}=0.981\pm 0.003$, which correspond to the normalised
stiffness ${\lambda}^{\text{est}}=95\pm 15$ [$s^{-1}$]. There was no
dependence between the obtained values, the sampling frequency and exposure
time, which proves that the ARMA fitting allows for the estimation of the trap
stiffness free of the studied effects.

Figure 4.6: MA(1) coefficient for different sampling frequencies and video
camera exposure times. The level $\theta_{1}=0$ is also drawn.

Given estimates $\phi^{\text{est}}_{1}$ and $\theta^{\text{est}}_{1}$ we can
inverse the AR and MA filters, effectively removing their influence from the
data.

The parameter $\phi_{1}$ is determined by the stiffness of the harmonic trap.
For a given trap with known stiffness, we can apply the filter
$L=1-\phi^{\text{est}}_{1}\mathcal{S}_{-1}$ to the data, i.e. analyse
$X_{k}-\phi^{\text{est}}_{1}X_{k-1}$. This series is an MA part of the ARMA
model which provides insight into distortions caused by the experimental
setup. If the covariance function of $X_{k}-\phi^{\text{est}}X_{k-1}$ has two
non-zero values (for $h=0$ and $h=1$), it is confirmed that some short-time
disturbance is present in the system. It can be mainly additive noise for the
positive estimate of $r_{X}(1)$ or mainly blur for negative value of
$r_{X}(1)$. We have not observed non-zero values of the covariance for larger
values of $h$, see Fig. 4.7. There we present the estimated covariance
function normalised. i.e. divided by the assessed value of the variance; such
memory function is often called correlation function.

The presence of $h>1$ non-zero values of the covariance would suggest the
presence of a blur between few adjacent frames or a composition of a blur and
an additive noise as in Eq. 4.4. The lack of this type of memory confirms that
ARMA(1,1) is a model sufficient to fully describe the physical dynamics and
distortions.

Figure 4.7:  Estimated correlation function of the series
$X_{k}-\phi^{\text{est}}_{1}X_{k-1}$, the used estimate for $\phi_{1}$ was
$\phi^{\text{est}}_{1}=0.9816$.

If we have calibrated the setup, i.e. we have determined the MA part of the
model, we can retrieve the original undisturbed series by applying the filter
$(1+\theta^{\text{est}}_{1}\mathcal{S}_{-1})^{-1}$. This operator has the form
of the infinite series
$\sum_{k=0}^{\infty}(-\theta^{\text{est}}_{1})^{k}\mathcal{S}_{-k}$, so for
practical purposes only the first few terms need to be used. The error that we
make using this approximation decays exponentially with the number of terms
taken along, so it is easy to obtain sufficient accuracy. If
$\theta^{\text{est}}_{1}\approx\theta_{1}$, the retrieved series
${X}^{\text{org}}_{k}\mathrel{\mathop{:}}=(1+\theta^{\text{est}}_{1}\mathcal{S}_{-1})^{-1}X_{k}$
satisfies the original AR equation

$\displaystyle{X}^{\text{org}}_{k}-\phi_{1}{X}^{\text{org}}_{k-1}=(1+\theta^{\text{est}}_{1}\mathcal{S}_{-1})^{-1}(X_{k}-\phi_{1}X_{k-1})$
$\displaystyle=(1+\theta^{\text{est}}_{1}\mathcal{S}_{-1})^{-1}(1+\theta_{1}\mathcal{S}_{-1})\eta_{k}\approx\eta_{k}.$
(4.94)

In the case of blur, the calculated $\eta_{k}$ is exactly the original
physical white noise; for the case of the additive white noise it has larger
variance than the original process (equal to
$\frac{D^{2}}{2\lambda}\left(\phi_{1}^{-2}-1\right)$) which can be accounted
for by a proper rescaling.

The memory function of the exemplary calculated ${X}^{\text{org}}$ is shown in
Fig. 4.8. We have taken the same time series as in Fig. 4.7 and using the
estimated value $\theta^{\text{est}}_{1}=0.1795$ we applied the first 20 terms
of the inverse operator $(1+\theta^{\text{est}}_{1}\mathcal{S}_{-1})^{-1}$.
The result is an AR(1) process, which is confirmed by Fig. 4.8. The first
estimated value of pacf is $0.918$, so is exactly equal to the previously
obtained estimate $\theta^{\text{est}}_{1}$.

Figure 4.8: Estimated pacf of the approximated series
$(1+\theta^{\text{est}}_{1}\mathcal{S}_{-1})^{-1}X_{k}$,
$\theta^{\text{est}}_{1}=0.1795$.

We can apply this procedure to the measured time-series if only we know the
value of MA parameter $\theta_{1}$ for a given calibration of the used
equipment, effectively removing its influence. This can be a powerful
technique if we calibrate the tweezers analysing their behaviour in the
controlled conditions, e.g. for spherical beads embedded in water, and later
remove the disturbances for data acquired in more complex systems such as
biological media.

For our data, the pacf for low frequencies fits to the model, see Fig 4.9. By
two dashed lines near zero we denote the level of an expected statistical
error. Values between these lines are statistically insignificant.

Figure 4.9: Partial autocorrelation estimated from the bead’s trajectory for
low frequency of sampling $10^{3}$ [fps].

However, the high frequency data have the pacf of different type, see Fig
4.10. It looks like a geometric series with a negative rate which indicates
the proper adjustment of the model.

Figure 4.10: Partial autocorrelation estimated from the bead’s trajectory for
high frequency of sampling $10^{4}$ [fps].

It appears that a proper modification explaining the data is addition of the
first-order Moving Average ( MA(1) ) part to the AR(1) model. The MA(1) part
represents a process for which the present value depends linearly on one past
value of the external noise (not the process itself), i.e.,
$\xi_{k}+\theta_{1}\xi_{k-1}$. The complete model, with both AR(1) and MA(1)
parts, is expressed explicitly as

$X_{k}=\phi_{1}X_{k-1}+\xi_{k}+\theta_{1}\xi_{k-1}.$ (4.95)

Unfortunately, from this form it is not clear what is the direct dependence
between $X_{k}$ and the past values of the process, because the value
$\xi_{k-1}$ is not orthogonal to these values. But, using the identity
$X_{k-1}=\phi_{1}X_{k-2}+\xi_{k-1}+\theta_{1}\xi_{k-2}$ we may write
$\xi_{k-1}=X_{k-1}-\phi_{1}X_{k-2}-\theta_{1}\xi_{k-2}$. Substituting this
equality into Eq. 4.95 we obtain

$X_{k}=(\phi_{1}+\theta_{1})X_{k-1}-\phi_{1}\theta_{1}X_{k-2}-\theta_{1}^{2}\xi_{k-2}.$
(4.96)

Now, we may repeat this procedure using the identity for $\xi_{k-2}$ and the
above formula. Continuing recursively we obtain

$X_{k}=(\phi_{1}+\theta_{1})\sum_{j=1}^{\infty}(-\theta_{1})^{j-1}X_{k-j}.$
(4.97)

This formula shows the explicit form of a dependence between $X_{k}$ and the
past values $X_{k-j}$, therefore justifying why the partial autocorrelation
has the form of a geometric series fitting the behaviour of the analysed data.
The parameters $\phi_{1}$ and $\theta_{1}$ can be estimated using least-
squares or maximum-likelihood techniques. The parameter $\theta_{1}$ is not
unique, there are two possible choices, see end of the last section. In our
context the one with $0<\theta_{1}<1$ has a proper interpretation and this one
was chosen. Performed fit returned $\theta^{\text{est}}_{1}=0.212\pm 0.005$
and $\phi^{\text{est}}_{1}=0.976\pm 0.001$, which perfectly corresponds to the
ARMA(1,1) process with the pacf drawn in Fig. 4.10. Given uncertainties of
$\phi^{\text{est}}_{1}$ and $\theta^{\text{est}}_{1}$ are standard deviations
estimated from a sample of ten trajectories with 60 000 observations in each
of them. The obtained sample of ten estimates of $\phi_{1}$ can be considered
Gaussian on standard level of significance 0.05 by Shapiro-Wilk, Anderson-
Darling and Pearson $\chi^{2}$ tests, the same is true for $\theta_{1}$. This
agrees with the well-known theory of long-trajectory asymptotics of the used
estimators (Brockwell and Davis, 2006). These standard deviations are only
slightly greater than the deviations for these estimators for such an
ARMA(1,1) process (which was checked by Monte Carlo simulation). It means that
measurement imperfections do not distort the precision of estimation.

Of course, there is no possibility to rule out that performed estimation is
precise, but biased. However, note that the obtained value
$\phi^{\text{est}}_{1}$ is consistent with the measurements of the stiffness
for lower frequencies.

The ARMA(1,1) can be naturally explained as an influence of the high-frequency
CCD camera. When the frequency of the sampling is high, the CCD matrix has no
time to fully refresh between subsequent photos. The remainder of the last
frame is still visible on the current one (the recording is smudgy and
blurred), which causes the centre of the mass position $X_{k}$ to include part
of the value $X_{k-1}$ from the last frame. Instead of $X_{k}$ we observe
$X^{\prime}_{k}=X_{k}+\theta_{1}X_{k-1}$. Parameter $\theta_{1}$ is exactly
the amount of intensity left from the last frame on the actual one,
$0<\theta_{1}<1$. In this situation

$\displaystyle X^{\prime}_{k}$ $\displaystyle=X_{k}+\theta_{1}X_{k-1}$
$\displaystyle=\phi_{1}(X_{k-1}+\theta_{1}X_{k-2})+\xi_{k}+\theta_{1}\xi_{k-1}$
$\displaystyle=\phi_{1}X^{\prime}_{k-1}+\xi_{k}+\theta_{1}\xi_{k-1},$ (4.98)

so $X^{\prime}_{k}$ is the ARMA(1,1) process. It is straightforward to obtain
the covariance function of this new process

$\displaystyle r_{X^{\prime}}(j)$
$\displaystyle=\operatorname*{{cov}}(X_{k}+\theta_{1}X_{k-1},X_{k+j}+\theta_{1}X_{k-1+j})$
$\displaystyle=(1+\theta_{1}^{2})r_{X}(j)+\theta_{1}\,r_{X}(j-1)+\theta_{1}\,r_{X}(j+1).$
(4.99)

As we see it is a sum of three geometric sequences. The form of the revised
power spectral density $s_{X}$ follows from Proposition 4.3.2

$s_{X}(\omega)=c\frac{1+\theta_{1}^{2}+2\theta_{1}\cos(\omega)}{1+\phi_{1}^{2}-2\phi_{1}\cos(\omega)}.$
(4.100)

This revised psd fits perfectly the data, see Fig. 4.11. Hence, the ARMA(1,1)
in the case of high-frequency recordings describes better the experimental
data than the commonly used AR(1) model.

Figure 4.11: Estimate of the psd (grey line) and the fitted ARMA(1,1) psd
(dashed black line).

The previous considerations prove usefulness of the ARMA(1,1) model, but they
do not imply that the observed effect is due to the influence of the CCD
camera. To justify our claim we provide further experimental evidence. During
all our previous measurement the exposure time for taking a frame was 45 [$\mu
s$]. The longer exposition is in most of the situations used only when
necessary (e.g. in cases of low illumination), since the recording is becoming
more blurred. So, if our conclusion is correct, an increase of the exposition
time would incorporate an MA(1) part even for low frequencies. This indeed is
true, as seen in Fig. 4.12; compare this result to Figs. 4.9 and 4.10.

Figure 4.12: Estimated pacf of trajectory taken with frequency of sampling
$10^{3}$ [fps] but exposure time increased to 895 [$\mu s$].

The obtained effect is identical to the MA(1) part (with coefficient
$\theta_{1}^{\text{est}}=0.172\pm 0.005$) present for high frequency and low
exposure time data. Analogical operation performed for larger sampling
frequencies increases the value of the MA(1) coefficients. Unfortunately, a
decrease of the exposition, which would likely diminish the MA(1) part for
higher frequencies, is hard to obtain from technical reasons, which leaves
described methodology as the only practical way to deal with this influence.

##  Chapter 5 Ergodicity and Fourier space

The relation between time averages and ensemble averages is still one of the
most important topics of statistical physics and this area of research is
under intense development. Its mathematical description, abstract ergodic
theory is a very wide subject (Walters, 1982; Cornfel et al., 1982). However,
in recent years, a new trend has emerged which concentrates on very practical
questions, concerning the behaviour of time-averaged statistics useful in
analysis of the real data. The most prominent example of a quantity used in
this approach is the time-averaged mean square displacement (Jeon et al.,
2011; Metzler et al., 2014).

In this chapter we use the Maruyma’s ergodic theorem (Maruyama, 1970) to study
the long-time behaviour of linear response systems, in particular modelled by
the classical and generalized Langevin equations. This theorem allow us to
formulate simple criteria which determine the ergodicity of the studied
processes and provides interesting insight into the physical origin of
ergodicity and non-ergodicity, which relates them to the structure of a heat
bath and the Hamiltonian of the system. It also allows to analyse Gaussian
models which are ergodic, but non-mixing which is an often overlooked
possibility.

Our most important result is a generalization of the Maruyma’s theorem
(Ślęzak, 2017), which determines the behaviour of the time-averages for
processes which are non-ergodic but stationary, and have profound physical
interpretation based on theory of harmonic oscillations, phonons and Langevin
equation. We demonstrate that, in the case of these processes, the disparity
between ensemble- and time-averages can be studied using non-linear
statistics, such as time-averaged characteristic function.

The obtained result is also interesting from the theoretical standpoint, as
its proof links mathematical concepts from the number theory, group theory and
second order stochastic processes, which is an unusual, yet insightful,
occurrence in the field of applied mathematics.

###  5.1 Basics of ergodic theory

Ergodic theory studies the behaviour of various averages of $f(X)$, where $f$
is a function of the whole trajectory $X$.

######  Definition 5.1.1.

For a given stochastic process $X$, if for a function $f$

$\mathbb{E}|f(X)|<\infty,$ (5.1)

we call this function an observable. If

$\mathbb{E}|f(X)|^{2}<\infty,$ (5.2)

we call $f$ a second-order observable.

Examples include observable of mean position $f(X)=X(t)$, mean square
displacement $f(X)=\big{(}X(t+\Delta)-X(t)\big{)}^{2}$, covariance
$f(X)=X(t+\Delta)X(t)$, and others. Take note, that for the above examples to
make sense, the measured averages must be time independent which will indeed
be the case. That means we consider processes which are stationary, or at
least have stationary increments (see Definition 2.2.2). Under such
assumptions one may take $t=0$ without any lose of generality.

If the process has a time-varying mean
$m(t):=\mathbb{E}[X(t)]\neq\text{const}$ we can always decompose it as a
random, zero-mean part and deterministic non-zero part. Because we will study
systems in which complex behaviour will be contained in the random part, we
will assume $m(t)=0$; the case $m(t)\neq 0$ would be a straightforward
generalisation.

For every process there exists an associated family of time-shift operators
$\mathcal{S}_{\tau}$ which describe the temporal evolution of the system, i.e.
$\mathcal{S}_{\tau}X=(X(t+\tau))_{t\in\mathbb{R}}$. The time-shift operator
was already used as a tool to analyse discrete-time ARMA processes, now we
consider its continuous-time counterpart.

With that in mind we can formulate the ergodic theorem. We will use general
form of this result, which gives deeper insight into behaviour of non-ergodic
processes.

######  Theorem 5.1.1 (Birkhoff’s ergodic theorem).

For a stationary stochastic process $X$, that is, for a process with the
measure-preserving shift operators $\mathcal{S}_{\tau}$, and any observable
$f$

$\lim_{T\to\infty}\frac{1}{T}\int_{0}^{T}\mathrm{d}\tau\
f\big{(}\mathcal{S}_{\tau}X\big{)}=\mathbb{E}[f(X)|\mathcal{C}],$ (5.3)

where by $\mathbb{E}[\boldsymbol{\cdot}|\mathcal{C}]$ we understand
conditional expected value under condition $\mathcal{C}$, with $\mathcal{C}$
being the $\sigma$-algebra of sets invariant under family of transformations
$\\{\mathcal{S}_{\tau}\\}_{\tau\in\mathbb{R}}$.

Essentially we calculate the expected value assuming that all time invariant
properties of $X$ are fixed. The physical interpretation of $\mathcal{C}$ is
that this is a set of constants of motion associated with $X$.

By the measure preserving transformation we mean that

$\operatorname*{P}(f(X)\in A)=\operatorname*{P}(f(\mathcal{S}_{\tau}X)\in A)$
(5.4)

for all measurable $A$. This essentially means that $X$ and time shifted
$\mathcal{S}_{\tau}X$ are statistically undistinguishable, which coincides
with the previously introduced notion of stationarity.

Therefore the Birkhoff’s ergodic theorem states that if only the process $X$
is stationary, the time-average of any observable converges to a random
variable $\mathbb{E}[f(X)|\mathcal{C}]$ which we can precisely determine if we
can identify all time invariant properties of $X$.

It is now clear that for the classical ergodic theorem to hold, the process
should be stationary and the condition $\mathcal{C}$ must be sufficiently
weak, so that $\mathbb{E}[f(X)|\mathcal{C}]=\mathbb{E}[f(X)]$; there should be
no significant time-invariant properties of $X$.

One immediate and important consequence of the Birkhoff’s theorem is that,
because of the so-called tower property
$\mathbb{E}\big{[}E[\boldsymbol{\cdot}|\mathcal{C}]\big{]}=\mathbb{E}[\boldsymbol{\cdot}]$,
for a stationary process the ensemble average of the time average of any
observable is equal to the ensemble average

$\mathbb{E}\left[\lim_{T\to\infty}\frac{1}{T}\int_{0}^{T}\mathrm{d}\tau\
f\big{(}\mathcal{S}_{\tau}X\big{)}\right]=\mathbb{E}[f(X)],$ (5.5)

in particular there is no possibility for a phenomena such as the weak
ergodicity breaking

$\mathbb{E}\left[\overline{\delta}{}^{2}_{X}(t)\right]\neq\mathbb{E}\left[\big{(}X(\tau+t)-X(\tau)\big{)}^{2}\right]$
(5.6)

studied e.g. in (Metzler et al., 2014).

Another concept closely related to the ergodicity is mixing. We say that the
system is mixing if for any observables $f,g$ such that
$\mathbb{E}|f(X)g(X)|<\infty$,

$\lim_{T\to\infty}\mathbb{E}[f(X)g(\mathcal{S}_{T}X)]=\mathbb{E}[f(X)]\mathbb{E}[g(X)].$
(5.7)

This is basically a statement that the process $X$ and the shifted process
$\mathcal{S}_{T}X$ are asymptotically independent; the dynamics of $X$ leaves
no persisting memory. It implies ergodicity and is often easier to study than
the ergodicity itself; however, there are examples of non-mixing ergodic
processes, even in the Gaussian case, which we will show in Section 5.4.

###  5.2 Classical Maruyama’s theorem

The main part of our considerations is true only for the class of Gaussian
processes. The precise definition 2.2.1 was stated in Chapter 2, but
intuitively speaking the sufficient and necessary condition for the process to
be Gaussian is that all $X(t)$ are Gaussian (we also admit the case when the
process is non-random as a degenerate Gaussian with variance 0) and they are
only linearly dependent (Janson, 1997). The presence of non-linear dynamics
excludes Gaussianity. This topic will be discussed in more detail in Chapter
6.

The limitation on the possible memory type has large consequences: the
Gaussian variables are fully described by their linear dependence structure
which is reflected in second moments. For stationary Gaussian process the
description is further simplified, because then
$m_{X}(t)=\mathbb{E}[X(t)]=\text{const}$ and $r_{X}(s,t)=r_{X}(t-s)$
(Proposition 2.27). Non-zero constant mean does not change any significant
properties of the process, so further on in this chapter we will consider only
zero-mean case. The above facts, and the linear structure of Gaussian
processes allow the mixing condition to be greatly reduced.

######  Proposition 5.2.1.

A stationary Gaussian process is mixing if and only if

$\lim_{T\to\infty}r_{X}(T)=0.$ (5.8)

This condition is often straightforward to check.

###### Proof.

If the process is mixing, taking $f(X)=g(X)=X(t)$ we obtain

$\lim_{T\to\infty}r_{X}(T)=\lim_{T\to\infty}\mathbb{E}[f(X)g(\mathcal{S}_{T}X)]=\mathbb{E}[f(X)]\mathbb{E}[g(X)]=0.$
(5.9)

To prove another implication we first show that

$\lim_{T\to\infty}\mathbb{E}[f(X)f(\mathcal{S}_{T}X)]=\big{(}\mathbb{E}[f(X)]\big{)}^{2}$
(5.10)

for any $f$ such that $\mathbb{E}|f(X)|^{2}<\infty$. It is sufficient to show
this property for the subclass of functions which are dense in the space of
observables. A natural choice is the space of finite-dimensional
characteristic functions, which contains the full information about the
distribution of $X$

$f(X)=\exp\left(\mathrm{i}\sum_{j=1}^{L}\theta_{j}X(t+t_{j})\right).$ (5.11)

The considered expected value is

$\displaystyle\mathbb{E}[f(X)f(\mathcal{S}_{T}X)]=\mathbb{E}\left[\exp\left(\mathrm{i}\sum_{j=1}^{L}\theta_{j}\big{(}X(t+t_{j})+X(t+t_{j}+T)\big{)}\right)\right]$
$\displaystyle=\exp\left(\sum_{i,j=1}^{L}\theta_{i}\theta_{j}\big{(}2r_{X}(t_{j}-t_{i})+r_{X}(t_{j}-t_{i}+T)+r_{X}(t_{j}-t_{i}+T)\big{)}\right).$
(5.12)

Taking the limit $T\to\infty$ we can use the fact that $r_{X}(t+T)\to 0$ and
get

$\lim_{T\to\infty}\mathbb{E}[f(X)f(\mathcal{S}_{T}X)]=\exp\left(2\sum_{i,j=1}^{L}\theta_{i}\theta_{j}r_{X}(t_{j}-t_{i})\right)=\big{(}\mathbb{E}[f(X)]\big{)}^{2}.$
(5.13)

Now, consider the observable $f+g$. Using the above result and mixing,

$\displaystyle\lim_{T\to\infty}\mathbb{E}[(f(X)+g(X))(f(\mathcal{S}_{T}X)+g(\mathcal{S}_{T}X))]=\big{(}\mathbb{E}[f(X)]+\mathbb{E}[g(X)]\big{)}^{2}$
$\displaystyle=\big{(}\mathbb{E}[f(X)]\big{)}^{2}+2\mathbb{E}[f(X)]\mathbb{E}[g(X)]+\big{(}\mathbb{E}[g(X)]\big{)}^{2}$
$\displaystyle=\lim_{T\to\infty}\big{(}\mathbb{E}[f(X)f(\mathcal{S}_{T}X)]+\mathbb{E}[f(X)g(\mathcal{S}_{T}X)]+\mathbb{E}[g(X)f(\mathcal{S}_{T}X)]+\mathbb{E}[g(X)g(\mathcal{S}_{T}X)]\big{)}$
$\displaystyle=\big{(}\mathbb{E}[f(X)]\big{)}^{2}+\lim_{T\to\infty}\big{(}\mathbb{E}[f(X)g(\mathcal{S}_{T}X)]+\mathbb{E}[g(X)f(\mathcal{S}_{T}X)]\big{)}+\big{(}\mathbb{E}[g(X)]\big{)}^{2}.$
(5.14)

If we compare the second and fourth line we get the equality

$\lim_{T\to\infty}\big{(}\mathbb{E}[f(X)g(\mathcal{S}_{T}X)]+\mathbb{E}[g(X)f(\mathcal{S}_{T}X)]\big{)}=2\mathbb{E}[f(X)]\mathbb{E}[g(X)].$
(5.15)

The two terms on the left are in fact equal. These expected values depend only
on the distribution of $X$. But the distribution of the process reflected in
time and process shifted in time is identical to the original one, so

$\mathbb{E}[g(X)f(\mathcal{S}_{T}X)]=\mathbb{E}[g(X)f(\mathcal{S}_{-T}X)]=\mathbb{E}[g(\mathcal{S}_{T}X)f(X)]=\mathbb{E}[f(X)g(\mathcal{S}_{T}X)].$
(5.16)

Finally

$\lim_{T\to\infty}\mathbb{E}[f(X)g(\mathcal{S}_{T}X)]=\mathbb{E}[f(X)]\mathbb{E}[g(X)].$
(5.17)

∎

The ergodicity itself can also be expressed in the language of the covariance
function. Initially, the notion of metric transitivity was used in this
context (Grenander, 1950; Dym and McKean, 1976). The most important theorem in
this field was proven by Gishiro Maryuama (Maruyama, 1970).

######  Theorem 5.2.1 (Maruyama).

A Gaussian stationary process $X$ is ergodic if and only if

$\lim_{T\to\infty}\frac{1}{T}\int_{0}^{T}\mathrm{d}\tau\ |r_{X}(\tau)|=0.$
(5.18)

The presence of the modulus $|r_{X}(\tau)|$ is crucially important, because it
excludes periodic oscillations of $r_{X}$. Generally this condition may also
be easy to check, as it is enough to know the asymptotic tail behaviour of the
covariance function. But, at the same time, this theorem does not give much
insight into the memory structure of non-ergodic Gaussian processes which are
studied in Section 5.5. In fact it can be viewed as a consequence of the more
general theorem shown there.

###  5.3 Harmonic processes

Consider elementary example of a motion in the harmonic potential, governed by
the equation

$\frac{\mathrm{d}}{\mathrm{d}t}{X}=-\omega_{0}^{2}X,\quad
X(0)=X_{0},\frac{\mathrm{d}}{\mathrm{d}t}{X}(0)=V_{0},$ (5.19)

with has the solution

$X(t)=X_{0}\cos(\omega_{0}t)+\frac{V_{0}}{\omega_{0}}\sin(\omega_{0}t).$
(5.20)

Under the assumption that at the beginning of the evolution, the system
interacted with the heat bath, $X_{0}$ and $V_{0}$ have Gibbs distribution
given by the density

$p(x_{0},v_{0})\propto\exp\left(-\omega_{0}^{2}\frac{x_{0}^{2}}{2k_{B}\mathcal{T}}\right)\exp\left(-\frac{v_{0}^{2}}{2k_{B}\mathcal{T}}\right).$
(5.21)

The resulting stochastic process has the covariance function

$r_{X}(t)=\frac{k_{B}\mathcal{T}}{\omega_{0}^{2}}\cos(\omega_{0}t),$ (5.22)

so it is a stationary Gaussian process. Its spectral representation is
therefore given by the measure

$\sigma_{X}(\mathrm{d}\omega)=\frac{k_{B}\mathcal{T}}{2\omega_{0}^{2}}\big{(}\delta(\mathrm{d}\omega-\omega_{0})+\delta(\mathrm{d}\omega+\omega_{0})\big{)}$
(5.23)

which is concentrated in 2 points, which we denote by two Dirac deltas. In a
natural way a question about ergodicity arises. Whereas the time average of
the observable of the position $f(X)=X(t)$

$\lim_{T\to\infty}\frac{1}{T}\int_{0}^{T}\mathrm{d}\tau\
X(\tau)=\lim_{T\to\infty}\frac{1}{T}\left(\frac{X_{0}}{\omega_{0}}\sin(\omega_{0}T)=\frac{V_{0}}{\omega_{0}^{2}}(\cos(\omega_{0}T)-1)\right)=0$
(5.24)

converges to the ensemble mean $0=\mathbb{E}[X(t)]$, the observable of the
mean square displacement does not, as

$\lim_{T\to\infty}\frac{1}{T}\int_{0}^{T}\mathrm{d}\tau\
\big{(}X(\tau+\Delta)-X(\tau)\big{)}^{2}=2\left(X_{0}^{2}+\frac{V_{0}^{2}}{\omega_{0}^{2}}\right)\left(\sin\left(\frac{\omega_{0}\Delta}{2}\right)\right)^{2}$
(5.25)

differs from the ensemble average

$\mathbb{E}\left[\big{(}X(t+\Delta)-X(t)\big{)}^{2}\right]=4\frac{k_{B}\mathcal{T}}{\omega_{0}^{2}}\left(\sin\left(\frac{\omega_{0}\Delta}{2}\right)\right)^{2}.$
(5.26)

From the point of view of Hamiltonian mechanics this lack of ergodicity is
expected; after the initial contact with a heat bath the system evolves as a
microcanonical ensemble and the trajectories are trapped on the surface of
constant energy, which prohibits ergodicity. Indeed, the term
$X_{0}^{2}+V_{0}^{2}/\omega_{0}^{2}$ on the left side of Eq. 5.25 is the total
energy of the system, generally random, but constant on the trajectories of
$X$; whereas $2k_{B}\mathcal{T}/\omega_{0}^{2}$ in 5.26 is the mean total
energy.

A similar reasoning applies to the more general process of the form

$X(t)=\sum_{k}A_{k}\mathrm{e}^{\mathrm{i}\omega_{k}t},$ (5.27)

where the sum can even be infinite if $A_{k}$ are independent, complex
Gaussian variables and $\sum_{k}\mathbb{E}|A_{k}|^{2}<\infty$. Random
functions of this class, called harmonic processes, are appearing e.g. in
phonon theory. Their covariance function and spectral measure are

$r_{X}(t)=\sum_{k}\mathbb{E}|A_{k}|^{2}\cos(\omega_{k}t),\quad\sigma_{X}(\mathrm{d}\omega)=\sum_{k}\frac{\mathbb{E}|A_{k}|^{2}}{2}\big{(}\delta(\mathrm{d}\omega-\omega_{k})+\delta(\mathrm{d}\omega+\omega_{k})\big{)}.$
(5.28)

If one calculates the ensemble- and time-average of the mean-square
displacement for such a process, the different nodes of oscillation prove to
be uncoupled in both time- and ensemble- average sense; the corresponding
formulas are sums of terms as in Eq. 5.26 or 5.25. Therefore, any process
within this class is stationary, but non-ergodic. The next section will show
that it is the only case of non-ergodic Gaussian stationary process.

###  5.4 Spectral form of Maruyama’s theorem

All the properties of a Gaussian process can be described interchangeably by
its covariance function or its spectral measure; this very specific property
of the Gaussian class is caused by its linear structure. For sure Maruyama’s
theorem can be expressed in the language of the spectral measure, and this
reformulation leads to a surprisingly elegant statement (Dym and McKean, 1976)

######  Theorem 5.4.1.

A stationary Gaussian process is ergodic if and only if its spectral measure
has no atoms.

Any measure can be decomposed as a sum of three distinct components: the
absolutely continuous, singular and discrete measures. For a stochastic
process the corresponding decomposition of the spectral measure
$\sigma=\sigma_{\text{ac}}+\sigma_{\text{s}}+\sigma_{\text{d}}$ causes also
the process itself to decompose into three independent components

$X(t)=X_{\text{ac}}(t)+X_{\text{s}}(t)+X_{\text{d}}(t),$ (5.29)

which is guaranteed by the harmonic representation 2.35.

  * i)

The component $X_{\text{ac}}$ is mixing. It has a power spectral density. The
Riemann-Lebesgue lemma shows that in this situation the covariance of
$X_{\text{ac}}$, and all memory functions decay at infinity, i.e. the values
of the process become asymptotically independent at long time scales.

  * ii)

The component $X_{\text{s}}$ is ergodic, but its memory structure may be
complex. Its covariance function does not necessarily decay to 0. It may
oscillate, but must be aperiodic and the high correlation events must become
more scarce as $t\to\infty$.

  * iii)

The last component $X_{\text{d}}$ is non-ergodic and is a Gaussian harmonic
process.

The case ii) is most mathematically challenging, because the set of spectral
measures, for which the corresponding covariance function decays, called
Rajchman measures, does not have a convenient description (Lyons, 1995). As a
demonstration let us consider an example using the best-known singular
measure: the Cantor measure.

The Cantor set is obtained by removing from the middle one-third the interval
$[0,1]$, then repeating this procedure for the two remaining intervals
$[0,1/3]$, $[1/3,1]$ and recursively applying this procedure infinitely many
times. The points which will remain are Cantor points. They can be more easily
characterised as points in $[0,1]$ which have no 1s in their ternary
representation, i.e. have the representation
$\sum_{k=1}^{\infty}d_{k}3^{-k},d_{k}\in\\{0,2\\}$.

The Cantor measure $\sigma_{C}$ is the uniform measure on the Cantor set. We
move it by -1/2 to the interval $[-1/2,1/2]$ so that the corresponding process
will be real-valued. Elementary calculation proves that the length of the
intervals removed during construction of the Cantor set is 1, therefore the
Cantor measure cannot have a density and must be singular. It is probably the
simplest to understand it as a discrete uniform distribution on the i.i.d.
series $D_{k}$, $\operatorname*{P}(D_{k}=-1)=\operatorname*{P}(D_{k}=1)=1/2$
mapped onto interval $[-1/2,1/2]$ using the formula
$Y=\sum_{k=1}^{\infty}D_{k}3^{-k}$. The process $C$ which has the Cantor
measure as a spectral measure has the covariance function

$r_{C}(t)=\int_{\mathbb{R}}\sigma_{C}(\mathrm{d}\omega)\
\mathrm{e}^{\mathrm{i}\omega
t}=\mathbb{E}\left[\mathrm{e}^{\mathrm{i}t\sum_{k=1}^{\infty}D_{k}3^{-k}}\right]=\prod_{k=1}^{\infty}\mathbb{E}\left[\mathrm{e}^{\mathrm{i}tD_{k}3^{-k}}\right]=\prod_{k=1}^{\infty}\cos\left(t3^{-k}\right).$
(5.30)

In contrast to the better-known classes of the covariance functions, $r_{C}$
has a specific property close to a self-similarity

$r_{C}(3t)=\cos(t)r_{C}(t),$ (5.31)

which also guarantees that $r_{C}$ does not decay to zero. The extremal points
of $r_{C}$ are located at $t_{k}=3^{k}\pi$,
$r_{C}(t_{k})=(-1)^{k+1}r_{C}(\pi)$, where it attains the values $\approx\pm
0.47$, see Fig. 5.1. It may not be clear that the function $r_{C}$ can be
easily calculated numerically, but the convergence is actually quite fast.

######  Proposition 5.4.1.

For the numerically approximated covariance function of the Cantor process

$\tilde{r}(t)\mathrel{\mathop{:}}=\prod_{k=1}^{N}\cos\big{(}3^{-k}t\big{)},$
(5.32)

the above function converges to the Cantor covariance function $r_{C}$ and for
$N\geq\log_{3}t$ it is bounded by

$(1+ct^{2}9^{-N})r_{C}(t)\leq\tilde{r}(t)\leq\mathrm{e}^{dt^{2}9^{-N}}r_{C}(t),\quad
c,d>0.$ (5.33)

This inequality holds for those values of $t$ for which $r_{C}(t)$ is
positive, and the reverse inequality holds for those $t$ when it is negative.

###### Proof.

For $N\geq\log_{3}t$ and any $j\in\mathbb{N}$ we have
$x=3^{-(N+j)}t\leq\pi/2$, $\cos\big{(}3^{-(N+j)}t\big{)}>0$ and

$1-x^{2}/2\leq\cos(x)\leq 1-4/\pi^{2}x^{2}.$ (5.34)

Expressing $\tilde{r}$ in terms of $r_{C}$

$\tilde{r}(t)=r_{C}(t)\prod_{j=1}^{\infty}\frac{1}{\cos\big{(}3^{-(N+j)}t\big{)}}$
(5.35)

and assuming that $r_{C}$ is positive for a given $t$,

$r_{C}(t)\prod_{j=1}^{\infty}\frac{1}{1-3^{-2(N+j)}t^{2}4/\pi^{2}}\leq\tilde{r}(t)\leq
r_{C}(t)\prod_{j=1}^{\infty}\frac{1}{1-3^{-2(N+j)}t^{2}/2}.$ (5.36)

We need to approximate the product of terms

$\frac{1}{1-a9^{-N}9^{-j}}=1+\frac{1}{a^{-1}9^{N}9^{j}-1}=1+9^{-N}\frac{1}{a^{-1}9^{j}-9^{-N}}\mathop{=\colon}1+p_{j}.$
(5.37)

To use the inequality valid for positive $p_{j}$’s

$1+\sum_{j=1}^{\infty}p_{j}\leq\prod_{j=1}^{\infty}(1+p_{j})\leq\exp\left(\sum_{j=1}^{\infty}p_{j}\right)$
(5.38)

we estimate the sum of $p_{j}$’s by

$\sum_{j=1}^{\infty}p_{j}\geq a9^{-N}\sum_{j=1}^{\infty}\frac{1}{9^{j}}\geq
a9^{-N}\frac{1}{8}$ (5.39)

from the bottom, and writing

$\sum_{j=1}^{\infty}p_{j}=9^{-N}\frac{1}{a^{-1}-9^{-N}}+9^{-N}\sum_{j=1}^{\infty}\frac{1}{9a^{-1}9^{j}-9^{-N}}\leq\frac{1}{9^{N}a^{-1}-1}+\frac{1}{9}\sum_{j=1}^{\infty}p_{j}$
(5.40)

we estimate

$\sum_{j=1}^{\infty}p_{j}\leq\frac{9}{8}\frac{1}{9^{N}a^{-1}-1}$ (5.41)

from the top. Substituting the proper $a$ we obtain the final result

$\left(1+\frac{t^{2}}{2\pi^{2}}9^{-N}\right)r_{C}(t)\leq\tilde{r}(t)\leq\exp\left(\frac{9}{8}\frac{1}{2t^{-2}-9^{-N}}9^{-N}\right)r_{C}(t).$
(5.42)

For the negative $r_{C}(t)$ the above inequality is reversed. ∎

Figure 5.1: Plot of the covariance function $r_{C}$ of the process with the
Cantor spectral measure.

The process $(C(t))_{t}$ exhibits recurring correlations that prohibit mixing;
however, events of high-dependence are becoming exponentially rare as the time
delay increases, which allows for the emergence of ergodicity.

The situation becomes more complicated with even slight generalisations of the
model. If instead we perform the recursive removal procedure such that at any
step the remaining intervals on the left and right have the length
one-$\eta$th of the previous one ($\eta$ being real number bigger than $2$),
the obtained singular measure and the process is non-mixing for natural $n$,
but it is mixing for any $\eta$ which is not a Pisot-Vijayaraghavan number
(Salem, 1963). The Pisot-Vijayaroghavan numbers are a closed countable set
which causes even infinitesimally small changes of $\eta$ to change the mixing
behaviour.

The complex ergodic behaviour complicates the analysis of the models with
singular spectral measures, but it is worth stressing that the erratic
behaviour of the covariance functions may be useful for describing
observations which could be otherwise accounted for as experimental errors. It
is also worth noting that the singular measures are gaining attention for
their relation to fractal dynamics and self-similarity (Samorodnitsky and
Taqqu, 1994; Mandelbrot, 2002; Embrechts and Mejima, 2002).

###  5.5 Generalised Maruyama’s theorem

Any real stationary Gaussian process can be written as

$X(t)=X_{\text{erg}}(t)+\sum_{k=1}^{N}R_{k}\cos(\Theta_{k}+\omega_{k}t)+X_{0},\quad\omega_{k}\neq
0,$ (5.43)

which follows from the ergodic decomposition 5.29, after taking the real part
of the harmonic process 5.27. The variables
$R_{k}\mathrel{\mathop{:}}=|A_{k}|$ are amplitudes of the spectral points at
frequencies $\omega_{k}$ and have Rayleigh distribution with densities

$p_{R_{k}}(x)=\frac{x^{2}}{\sigma_{k}^{2}}\exp\left(-\frac{x^{2}}{2\sigma_{k}^{2}}\right)$
(5.44)

and scale parameters $\sigma_{k}^{2}=\mathbb{E}|A_{k}|^{2}$. Variables
$\Theta_{k}$ are phases of $A_{k}$ and have uniform distribution on
$[0,2\pi)$, a consequence of rotational invariance of i.i.d. Gaussian vectors.
In full generality the number of spectral points may be infinite $N=\infty$,
however in this case the process $X$ may exhibit a complicated aperiodic
behaviour; as it is not very important for the applicational purposes, most of
our subsequent result will be true only in the case $N<\infty$.

The decomposition into non-ergodic and ergodic components yields a useful and
straightforward description of the statistical properties of Gaussian
processes. It is made possible by the full characterisation of the invariant
sets of this dynamical system.

######  Theorem 5.5.1.

For any stationary Gaussian process $X$ with $N<\infty$ spectral points at
rationally incommensurable frequencies $\\{\omega_{k}\\}_{k=1}^{N}$, the
family of invariant sets $\mathcal{C}$ is the $\sigma$-algebra
$\sigma(X_{0},\\{R_{k}\\}_{k})$

###### Proof.

We will combine methods for trigonometric series presented in (Kimme, 1965)
and Gaussian ergodic theorem from Section 3.9 of (Dym and McKean, 1976).

We will use the representation

$X(t)=X_{\text{erg}}(t)+\sum_{k=1}^{N}R_{k}\cos(\Theta_{k}+\omega_{k}t)+X_{0},\quad\omega_{k}\neq
0,$ (5.45)

where the $\omega_{k}$ are distinct. The process $X_{\text{erg}}(t)$ has no
spectral points. The full distribution of $X$ is generated by the values
$X_{\text{erg}}(t_{l}),R_{j}\cos(\Theta_{j}+\omega_{j}t_{j}),X_{0}$. Therefore
it is sufficient to study the time-average distribution of the sum

$\sum_{l=1}^{L}\theta_{l}X_{\text{erg}}(t+t_{l})+\sum_{j=1}^{N}\lambda_{j}R_{j}\cos(\Theta_{j}+\omega_{j}t_{j}+\omega_{j}t)+\lambda_{0}X_{0}.$
(5.46)

For brevity we denote
$\tilde{\Theta}_{j}\mathrel{\mathop{:}}=\Theta_{j}+\omega_{j}t_{j}$. As the
distribution is uniquely determined by the corresponding characteristic
function, we will compare the time-averaged characteristic function and the
ensemble-averaged one given the condition $\sigma(X_{0},\\{R_{k}\\})$. Their
equality will prove the theorem.

First we calculate the conditional ensemble-averaged characteristic function.
The variables $\tilde{\Theta}_{j}$ have the same distribution as $\Theta_{j}$
modulo $2\pi$, because they are independent from each other, $R_{j}$’s,
$X_{0}$, $X_{\text{erg}}$, and have marginal uniform distribution. We get

$\displaystyle\
\mathbb{E}\left[\mathrm{e}^{\mathrm{i}\left(\sum_{l=1}^{L}\theta_{l}X_{\text{erg}}(t+t_{l})+\sum_{j=1}^{N}\lambda_{j}R_{j}\cos(\tilde{\Theta}_{j}+\omega_{j}t)+\lambda_{0}X_{0}\right)}|X_{0},\\{R_{k}\\}\right]$
$\displaystyle=\mathbb{E}\left[\mathrm{e}^{\mathrm{i}\left(\sum_{l=1}^{L}\theta_{l}X_{\text{erg}}(t_{l})\right)}\right]\prod_{j=1}^{N}\mathbb{E}\left[\mathrm{e}^{\mathrm{i}R_{j}\cos(\tilde{\Theta}_{j})}|R_{j}\right]\mathrm{e}^{\mathrm{i}\lambda_{0}X_{0}}$
$\displaystyle=\phi_{\theta_{1},\ldots,\theta_{L}}\prod_{j=1}^{N}J_{0}(R_{j}\lambda_{j})\mathrm{e}^{\mathrm{i}\lambda_{0}X_{0}},\quad\phi_{\theta_{1},\ldots,\theta_{L}}\mathrel{\mathop{:}}=\mathbb{E}\left[\mathrm{e}^{\mathrm{i}\left(\sum_{l=1}^{L}\theta_{l}X_{\text{erg}}(t_{l})\right)}\right],$
(5.47)

where $J_{0}$ are Bessel functions of the first kind and order $0$; they stem
from the formula

$\frac{1}{2\pi}\int_{-\pi}^{\pi}\mathrm{d}x\ \mathrm{e}^{\mathrm{i}\lambda
R\cos(x)}=J_{0}(\lambda R).$ (5.48)

The rest of the proof will be the calculation of the time-average.

We denote

$\Phi(t)\mathrel{\mathop{:}}=\mathrm{e}^{\mathrm{i}\left(\sum_{l=1}^{L}\theta_{l}X_{\text{erg}}(t+t_{l})\right)}$
(5.49)

which appears in the integral used during calculations of the time-average

$I_{T}\mathrel{\mathop{:}}=\int_{0}^{T}\mathrm{d}\tau\
\Phi(t)\mathrm{e}^{\mathrm{i}\left(\sum_{j=1}^{L}\lambda_{j}R_{j}\cos(\tilde{\Theta}_{j}+\omega_{j}t)\right)}\mathrm{e}^{\mathrm{i}\lambda_{0}X_{0}}.$
(5.50)

The factor $\mathrm{e}^{\mathrm{i}\lambda_{0}X_{0}}$ already agrees with the
conditional average, so we will assume $X_{0}=0$ later on for brevity.

Next we expand each exponent of cosine using the Jacobi-Anger identity

$\mathrm{e}^{\mathrm{i}z\cos(w)}=\sum_{m=-\infty}^{\infty}\mathrm{i}^{m}\mathrm{e}^{\mathrm{i}mw}J_{m}(z),$
(5.51)

obtaining

$\displaystyle I_{T}$ $\displaystyle=\int_{0}^{T}\mathrm{d}\tau\
\Phi(\tau)\prod_{j=1}^{N}\sum_{m=-\infty}^{\infty}\mathrm{i}^{m}\mathrm{e}^{\mathrm{i}m(\tilde{\Theta}_{j}+\omega_{j}\tau)}J_{m}(\lambda_{j}R_{j})$
$\displaystyle=\int_{0}^{T}\mathrm{d}\tau\ \Phi(\tau)\sum_{S\in
M_{N}}\prod_{m_{j}\in
S}\mathrm{e}^{\mathrm{i}m_{j}(\tilde{\Theta}_{j}+\omega_{j}\tau)}\mathrm{i}^{m_{j}}J_{m_{j}}(\lambda_{j}R_{j})$
$\displaystyle=\int_{0}^{T}\mathrm{d}\tau\ \Phi(\tau)\sum_{S\in
M_{N}}\exp\left(\mathrm{i}\sum_{m_{j}\in
S}m_{j}(\tilde{\Theta}_{j}+\omega_{j}\tau)\right)\prod_{m_{j}\in
S}\mathrm{i}^{m_{j}}J_{m_{j}}(\lambda_{j}R_{j})$ $\displaystyle=\sum_{S\in
M_{N}}\exp\left(\mathrm{i}\sum_{m_{j}\in
S}m_{j}\tilde{\Theta}_{j}\right)\prod_{m_{j}\in
S}\mathrm{i}^{m_{j}}J_{m_{j}}(\lambda_{j}R_{j})\int_{0}^{T}\mathrm{d}\tau\
\Phi(\tau)\exp\left(\mathrm{i}\tau\sum_{m_{j}\in S}m_{j}\omega_{j}\right)$
(5.52)

where $M_{N}$ is the family of all $N$-element subsets of integers $S$.
Exchanging the order of infinite sum and integral is possible because the
integrated function is bounded by $1$, which also justifies commuting the sum
and the limit in the next step. We shall denote
$\Omega_{S}\mathrel{\mathop{:}}=\sum_{m_{j}\in S}m_{j}\omega_{j}$ and check
that

$\lim_{T\to\infty}\frac{1}{T}\int_{0}^{T}\mathrm{d}\tau\
\Phi(\tau)\exp\left(\mathrm{i}\tau\Omega_{S}\right)=0,\text{ if
}\Omega_{S}\neq 0$ (5.53)

The proof of this statement is given in the lemma below. Let us conclude the
whole proof. Because $\\{\omega_{k}\\}$ are rationally incommensurable, the
equality $\Omega_{S}=\sum_{m_{j}\in S}m_{j}\omega_{j}=0$ holds for integer
$m_{j}$ only when $m_{1}=m_{2}=\ldots=m_{N}=0$. In the sum $\sum_{S\in M_{N}}$
only one element $S=\\{0,0,\ldots,0\\}$ remains and

$\lim_{T\to\infty}\frac{1}{T}I_{T}=\prod_{j=1}^{N}J_{0}(\lambda_{j}R_{j})\lim_{T\to\infty}\frac{1}{T}\int_{0}^{T}\mathrm{d}\tau\
\Phi(\tau)=\prod_{j=1}^{N}J_{0}(\lambda_{j}R_{j})\phi_{\theta_{1},\ldots,\theta_{L}},$
(5.54)

where the last equality holds due to the ergodicity of $X_{\text{erg}}$. This
is the desired conditional mean. ∎

######  Lemma 5.5.1.

$\lim_{T\to\infty}\frac{1}{T}\int_{0}^{T}\mathrm{d}\tau\
\Phi(\tau)\exp\left(\mathrm{i}\tau\Omega\right)=0,$ (5.55)

for $\Omega\neq 0$ and $\Phi$ given by 5.49.

###### Proof.

First note that $\Phi$ is a strictly stationary random process. Take
$U\overset{d}{=}\mathcal{U}([0,2\pi))$ independent from $\Phi$. The process
$t\mapsto\Phi(t)\exp(\mathrm{i}U+\mathrm{i}\Omega t)$ is also stationary and
has a finite first moment equal to 1. Therefore, the Birkhoff’s ergodic
theorem guarantees the time-average exists almost surely and equals

$\lim_{T\to\infty}\frac{1}{T}\int_{0}^{T}\mathrm{d}\tau\
\Phi(\tau)\exp\left(\mathrm{i}U+\mathrm{i}\tau\Omega\right)=\exp(\mathrm{i}U)\lim_{T\to\infty}\frac{1}{T}\int_{0}^{T}\mathrm{d}\tau\
\Phi(\tau)\exp\left(\mathrm{i}\tau\Omega\right)=\widehat{X}(\Omega).$ (5.56)

Thus, the limit 5.55 also exists almost surely and equals to a random variable
$\exp(-\mathrm{i}U)\widehat{X}(\Omega)$ (for more details on
$\widehat{X}(\Omega)$ see (Doob, 1990, Chapter XI.2, page 516)). We will prove
that it is 0.

Let $r$ be the covariance function of $X_{\text{erg}}$ and $\sigma$ its
continuous spectral measure. We will study
$\mathbb{E}|\boldsymbol{\cdot}|^{2}$ of the above time-average and show its
mean-square convergence to 0, which suffices to prove also the almost sure
convergence to the same limit:

$\displaystyle\
\frac{1}{T^{2}}\int_{0}^{T}\mathrm{d}\tau_{1}\int_{0}^{T}\mathrm{d}\tau_{2}\mathbb{E}\left[\exp\left(\mathrm{i}\sum_{j=1}^{L}\theta_{j}\big{(}X_{\text{erg}}(t_{j}+\tau_{1})-X_{\text{erg}}(t_{j}+\tau_{2})\big{)}\right)\right]\mathrm{e}^{\mathrm{i}\Omega(\tau_{1}-\tau_{2})}$
$\displaystyle=\frac{C}{T^{2}}\int_{0}^{T}\mathrm{d}\tau_{1}\int_{0}^{T}\mathrm{d}\tau_{2}\exp\left(\sum_{j,k=1}^{L}\theta_{j}\theta_{k}r(t_{k}-t_{j}+\tau_{1}-\tau_{2})\right)\mathrm{e}^{\mathrm{i}\Omega(\tau_{1}-\tau_{2})}$
$\displaystyle=\frac{C}{T^{2}}\int_{0}^{T}\mathrm{d}\tau_{1}\int_{0}^{T}\mathrm{d}\tau_{2}\exp\left(\int_{\mathbb{R}}\mathrm{d}\tilde{\sigma}(\omega)\mathrm{e}^{\mathrm{i}\omega(\tau_{1}-\tau_{2})}\right)\mathrm{e}^{\mathrm{i}\Omega(\tau_{1}-\tau_{2})},$
(5.57)

where we denoted by $C$ the factor before the integral and by $\tilde{\sigma}$
the modified spectral measure; it is just multiplied by a continuous function.

$C\mathrel{\mathop{:}}=\exp\left(2\sum_{j,k=1}^{L}\theta_{i}\theta_{j}r(t_{k}-t_{j})\right),\quad\mathrm{d}\tilde{\sigma}(\omega)\mathrel{\mathop{:}}=\mathrm{d}\sigma(\omega)\left|\sum_{j=1}^{L}\theta_{j}\mathrm{e}^{\mathrm{i}\omega
t_{j}}\right|^{2}$ (5.58)

Next, we expand external $\exp(\boldsymbol{\cdot})$ into the Taylor series,
obtaining

$\displaystyle\
\frac{C}{T^{2}}\int_{0}^{T}\mathrm{d}\tau_{1}\int_{0}^{T}\mathrm{d}\tau_{2}\left(1+\sum_{n=1}^{\infty}\frac{1}{n!}\left(\int_{\mathbb{R}}\mathrm{d}\tilde{\sigma}(\omega)\mathrm{e}^{\mathrm{i}\omega(\tau_{1}-\tau_{2})}\right)^{n}\right)\mathrm{e}^{\mathrm{i}\Omega(\tau_{1}-\tau_{2})}$
$\displaystyle=\frac{C}{T^{2}}\int_{0}^{T}\mathrm{d}\tau_{1}\int_{0}^{T}\mathrm{d}\tau_{2}\left(1+\sum_{n=1}^{\infty}\frac{1}{n!}\int_{\mathbb{R}}\mathrm{d}\tilde{\sigma}^{*n}(\omega)\mathrm{e}^{\mathrm{i}\omega(\tau_{1}-\tau_{2})}\right)\mathrm{e}^{\mathrm{i}\Omega(\tau_{1}-\tau_{2})},$
(5.59)

where $\tilde{\sigma}^{*n}$ is the $n$-fold convolution power of
$\tilde{\sigma}$. This Taylor series is uniformly convergent. We commute the
limit $T\to\infty$ with the sum and calculate the integrals; for the term
$n=1$ we have

$\frac{1}{T^{2}}\int_{0}^{T}\mathrm{d}\tau_{1}\int_{0}^{T}\mathrm{d}\tau_{2}\mathrm{e}^{\mathrm{i}\Omega(\tau_{1}-\tau_{2})}=\frac{1}{T^{2}}\frac{1}{\Omega^{2}}\left|\mathrm{e}^{\mathrm{i}\Omega
T}-1\right|^{2}\xrightarrow{T\to\infty}0,$ (5.60)

where the assumption $\Omega\neq 0$ is crucial. For any other term

$\displaystyle\
\lim_{T\to\infty}\frac{1}{T^{2}}\int_{0}^{T}\mathrm{d}\tau_{1}\int_{0}^{T}\mathrm{d}\tau_{2}\int_{\mathbb{R}}\mathrm{d}\tilde{\sigma}^{*n}(\omega)\mathrm{e}^{\mathrm{i}(\omega+\Omega)(\tau_{1}-\tau_{2})}$
$\displaystyle=\lim_{T\to\infty}\frac{1}{T^{2}}\int_{\mathbb{R}}\mathrm{d}\tilde{\sigma}^{*n}(\omega-\Omega)\int_{0}^{T}\mathrm{d}\tau_{1}\int_{0}^{T}\mathrm{d}\tau_{2}\mathrm{e}^{\mathrm{i}\omega(\tau_{1}-\tau_{2})}$
$\displaystyle=2\lim_{T\to\infty}\int_{\mathbb{R}}\mathrm{d}\tilde{\sigma}^{*n}(\omega-\Omega)\frac{1-\cos(\omega
T)}{(\omega T)^{2}}.$ (5.61)

In the last line one recognises the functional which returns the jump of the
measure $\tilde{\sigma}^{*n}$ at the point $\Omega$. But, the measure $\sigma$
is continuous and $\tilde{\sigma}^{*n}$ is also continuous; the result is
$\tilde{\sigma}^{*n}(\\{\Omega\\})=0$. ∎

The assumption that spectral points at $\\{\omega_{k}\\}_{k=1}^{N}$ are
rationally incommensurable means that they cannot be represented as
$\omega_{k}=q_{k}\alpha$ for any rational $q_{k}$’s. It is trivially fulfilled
in most real physical systems, in which $\omega_{k}$ are self-frequencies of
the harmonic oscillators and depend on the complex set of the system’s
parameters.

In such a case, the aforementioned theorem guarantees that for any observable
$f$, the time average converges to

$\lim_{T\to\infty}\frac{1}{T}\int_{0}^{T}\mathrm{d}\tau\
f(\mathcal{S}_{\tau}X)=\mathbb{E}\big{[}f(X)|,X_{0},\\{R_{k}\\}_{k}\big{]},$
(5.62)

i.e. to the ensemble average calculated under the condition that the
amplitudes of the spectral points and the constant term $X_{0}$ are fixed. For
an observable which depends on the one time moment of $X$ only,
$f(X)=f(X(t))$, the above formula simplifies to the explicit integral

$\displaystyle\lim_{T\to\infty}\frac{1}{T}\int_{0}^{T}\mathrm{d}\tau\
f(X(\tau))$
$\displaystyle=\int_{0}^{2\pi}\mathrm{d}\theta_{1}\int_{0}^{2\pi}\mathrm{d}\theta_{2}\ldots\int_{\mathbb{R}}\mathrm{d}x\frac{1}{\sqrt{2\pi}c}\mathrm{e}^{-\frac{x^{2}}{2c^{2}}}f\big{(}x+\sum_{k}R_{k}\cos(\theta_{k})+X_{0}\big{)},$
(5.63)

which depends only on the variance of the ergodic component
$c^{2}=\mathbb{E}\big{[}X_{\text{erg}}(t)^{2}\big{]}$, $R_{k}$’s and $X_{0}$,
which are random, but fixed for each trajectory.

In particular, using the cumulative distribution function method, one can
calculate the non-ergodic time-averaged probability density
$\overline{p}_{X_{k}}$ of any given discrete spectral component with amplitude
$R_{k}=R$

$\overline{p}_{X_{k}}(x)=p_{X_{k}}(x|R)=\frac{\big{(}\pi
R\big{)}^{-1}}{\sqrt{1-\big{(}x/R\big{)}^{2}}},\quad-R\leq x\leq R.$ (5.64)

This quantity should be observed, if one uses the time-average estimation of
the probability density, e.g. the kernel density estimators or the histogram
(Silverman, 1986).

This is also an interesting example of a singularity caused by the non-
ergodicity. The flat extrema of the cosine function are responsible for the
probability density divergence of type $x^{-1/2}$ at points $-R$ and $+R$.
However, this unusual behaviour is not easy to directly observe as typical
data contain an ergodic component, e.g. some kind of noise from the
experimental setup. In such a case, the observed empirical probability
distribution $\overline{p}$ is a convolution of the ergodic component’s
Gaussian density with a given stationary variance $\sigma^{2}$ and some number
of densities of type 5.64, see Figure 5.2. The singular concentration of the
probability mass around $-R$ and $+R$ distorts the tails of the distribution
in some specific way: it thickens them by moving the original distribution,
but thins them through dividing by a square root factor. More precisely, we
recognise that the convolution of the densities is an example of a Laplace
transform and use Abelian theorem (Postnikov, 1980) to obtain the asymptotic
behaviour of the left tail

$\displaystyle\overline{p}(-x)$ $\displaystyle=\frac{1}{\sqrt{2\pi^{3}}\sigma
R}\int_{-R}^{R}\mathrm{d}y\frac{1}{\sqrt{1-\big{(}x/R\big{)}^{2}}}\mathrm{e}^{-\frac{(y+x)^{2}}{2\sigma^{2}}}$
$\displaystyle=\mathrm{e}^{-\frac{(x-R)^{2}}{2\sigma^{2}}}\frac{1}{\sqrt{2\pi^{3}R}\sigma}\int_{0}^{2R}\mathrm{d}y\frac{1}{\sqrt{y}\sqrt{2-y/R}}\mathrm{e}^{-\frac{y^{2}}{2\sigma^{2}}}\mathrm{e}^{-\frac{y(x-R)}{\sigma^{2}}}$
$\displaystyle\sim\frac{1}{4\pi\sqrt{1-x/R}}\mathrm{e}^{-\frac{(x-R)^{2}}{2\sigma^{2}}},\quad
x\to\infty;$ (5.65)

symmetrically for the right tail. This result does not change significantly
for a finite number of spectral points, as it depends only on the presence of
singularities of the convolved densities.

If there is only one discrete spectral component, that is
$X(t)=X_{\text{erg}}(t)+R\cos(\Theta+\omega
t),\mathbb{E}[X_{\text{erg}}^{2}]=\sigma^{2}$, the time-averaged pdf can be
expressed as a series of Bessel functions,

$\displaystyle\overline{p}_{X}(x)$
$\displaystyle=\frac{1}{\sqrt{2\pi\sigma^{2}}}\frac{1}{2\pi}\int_{0}^{2\pi}\mathrm{d}\theta\
\mathrm{e}^{-\frac{(x-R\cos(\theta))^{2}}{2\sigma^{2}}}=\frac{1}{\sqrt{2\pi\sigma^{2}}}\mathrm{e}^{-\frac{x^{2}}{2\sigma^{2}}}\frac{1}{2\pi}\int_{0}^{2\pi}\mathrm{d}\theta\
\mathrm{e}^{\frac{xR\cos(\theta)}{\sigma^{2}}}\mathrm{e}^{-\frac{R^{2}\cos(\theta)^{2}}{2\sigma^{2}}}$
$\displaystyle=p_{X_{\text{erg}}}(x)\frac{1}{2\pi}\int_{0}^{2\pi}\mathrm{d}\theta\
\mathrm{e}^{\frac{xR\cos(\theta)}{\sigma^{2}}}\sum_{k=0}^{\infty}\frac{(-1)^{k}}{k!2^{k}}\frac{R^{2k}}{\sigma^{2k}}\cos(\theta)^{2k}$
$\displaystyle=p_{X_{\text{erg}}}(x)\sum_{k=0}^{\infty}\frac{(-1)^{k}}{k!2^{k}}\frac{\mathrm{d}^{2k}}{\mathrm{d}x^{2k}}\frac{1}{2\pi}\int_{0}^{2\pi}\mathrm{d}\theta\
\mathrm{e}^{\frac{xR\cos(\theta}{\sigma^{2}}}=p_{X_{\text{erg}}}(x)\sum_{k=0}^{\infty}\frac{(-1)^{k}}{k!2^{k}}\frac{\mathrm{d}^{2k}}{\mathrm{d}x^{2k}}I_{0}\left(\frac{xR}{\sigma^{2}}\right).$
(5.66)

By $I_{0}$ we denoted the modified Bessel function of the first kind and order
$0$; the derivatives of this function can be expressed using linear
combinations of $I_{n}$ if needed.

Both asymptotic formula 5.5 and series 5.5 differ considerably from the normal
distribution; for most realisations the amplitudes $R_{k}$ are large enough to
strongly affect the time-averaged probability density. This behaviour is
presented in Fig. 5.2. For nearly all realisations statistical tests also show
significant non-Gaussianity (e.g. Shapiro and Kolmogorov-Smirnov tests).
However, for small amplitudes of the $R_{k}$’s this effect could be less
noticeable.

Figure 5.2: Time-average kernel density estimation and ergodic density for a
Gaussian process with one discrete spectral component,
$\mathbb{E}\big{[}R_{k}^{2}\big{]}=\mathbb{E}\big{[}X_{\text{erg}}(t)^{2}\big{]}=1$.

The fast decay of the function $\exp(-x^{2})$ may complicate the analysis of
non-ergodicity through a density estimation. A more convenient method is to
use the time-averaged characteristic function. It is a time-average of the
observable $f(X)=\exp(\mathrm{i}\theta X(t))$, which for a non-ergodic
component equals

$\lim_{T\to\infty}\frac{1}{T}\int_{0}^{T}\mathrm{d}\tau\
\mathrm{e}^{\mathrm{i}\theta
X(\tau)}=\frac{1}{2\pi}\int_{-\pi}^{\pi}\mathrm{d}x\
\mathrm{e}^{\mathrm{i}\theta R\cos(x)}=J_{0}(R\theta),$ (5.67)

where $J_{0}$ is a Bessel function of the first kind and order $0$. Therefore
the time-averaged characteristic function $\overline{\phi}$ of any stationary
Gaussian process with incommensurable frequencies has form

$\overline{\phi}(\theta)=\mathrm{e}^{-c^{2}\theta^{2}/2}\prod_{k}J_{0}\big{(}R_{k}\theta\big{)}.$
(5.68)

Unlike for the ergodic case, this function has zeros determined by the zeros
of the function $J_{0}$, which may be approximated numerically, the first
being $x_{1}\approx 2.405$, the second $x_{2}\approx 5.520$. Therefore the
location of zeros for the time-averaged characteristic function may serve to
preliminary estimate the number and values of the amplitudes $R_{k}$. A more
precise estimation requires least-squares fitting with the function in Eq.
5.68. Exemplary results are shown in Fig. 5.3 which prove that it is a viable
method of estimation for a small number of spectral points.

The estimators of $R_{k}$ have the tendency to return some undershoot values
which cause negative bias, especially for lower lengths of the trajectories,
but are generally reliable. The results depend on the sample length $N$, but
do not depend on the sampling time $\Delta t$ as long as $\pi\Delta t$ is
incommensurable with $\\{\omega_{k}\\}_{k}$. If this is not true, the measured
time series has a periodic component and the proper values of the time-
averaged observables are obtained in the infill asymptotics, i.e. $\Delta t\to
0$, and $\Delta tN\to\infty$. Such requirements guarantee that the calculated
mean converges to the time integral $1/T\int_{0}^{T}\mathrm{d}\tau$ and
$T\to\infty$.

Figure 5.3: Estimation of $c=\mathbb{E}[X_{\text{erg}}(t)^{2}],R_{1}$ and
$R_{2}$ for different sample sizes $N$, using a least-square fit of the time-
averaged characteristic function for the process
$X_{\text{erg.}}(t)+3\cos(t)+2\cos(\sqrt{2}t)$.

Let us consider such models in which it may happen that $\\{\omega\\}_{k}$ are
commensurable. This situation may appear practically when $\\{\omega_{k}\\}$
by coincidence or due to some symmetry are close to being commensurable, that
is they can be expressed as $\omega_{k}=\alpha p_{k}/q_{k}+\epsilon_{k}$ where
$p_{k},q_{k}$ are small natural numbers and $\epsilon_{k}$ is small compared
to $1/T$.

For commensurable $\\{\omega_{k}\\}$ the harmonic process is in fact periodic
and the length of its period is proportional to the lowest common denominator
of $\omega_{k}$’s. We show an example in Fig. 5.4, where the empirical
characteristic function of the process
$R_{1}\cos(\Theta_{1}\omega_{1}t){+R_{2}\cos(\Theta_{2}+\omega_{2}t)}$ is
presented. For simplicity we fixed $\omega_{1}=1$ and $R_{1}=1,R_{2}=2$, as
the $R_{k}$’s are constants of motion in this case. Because of the
periodicity, the calculated time-average depends on the random initial phases
$\Theta_{1},\Theta_{2}$.

Figure 5.4: Time-averaged characteristic function of the process
$2\cos(\Theta_{1}+t)+2\cos(\Theta_{2}+\omega_{2}t)$ for different $\omega_{2}$
and 20 realisations of $\Theta_{1},\Theta_{2}$. The black line, corresponding
to the simulation for $\omega_{2}=\sqrt{2}$ up to numerical accuracy
$\epsilon\approx 10^{-16}$, agrees perfectly with ergodic average
$J_{0}(t)J_{0}(2t)$, see Eq. 5.68.

The exact dependence on both random amplitudes and angles is actually known;
the time-averaged characteristic function of any stationary Gaussian process
is given by the formula (see proof of the theorem below)

$\overline{\phi}(\theta)=\mathrm{e}^{-c^{2}\theta^{2}/2}\sum_{S\in
G}\exp\left(\mathrm{i}\sum_{m_{j}\in S}m_{j}\Theta_{j}\right)\prod_{m_{j}\in
S}\mathrm{i}^{m_{j}}J_{m_{j}}(\theta R_{j}),$ (5.69)

where $G$ is the family of sets of $m_{j}$’s for which
$\sum_{j}m_{j}\omega_{j}=0$. The result depends on linear combinations of the
random phases $\sum_{j}m_{j}\Theta_{j}\ (\mathrm{mod}\ 2\pi)$. This is true
not only for the time-averaged characteristic function, but any observable,
which is stated in the following theorem.

######  Theorem 5.5.2.

For any stationary Gaussian process $X$ with $N<\infty$ spectral points
$\\{\omega_{k}\\}_{k=1}^{N}$ (some of which may be commensurable), the family
of invariant sets $\mathcal{C}$ is the $\sigma$-algebra
$\sigma(X_{0},\\{R_{k}\\}_{k},\mathcal{M})$, where $X_{0}$ is the constant
term, $R_{k}$ are the amplitudes of the atoms of the spectral measure and
$\mathcal{M}$ is a family

$\mathcal{M}=\left\\{\sum_{m_{j}}m_{j}\Theta_{j}\ (\mathrm{mod}\
2\pi)\colon\sum_{j}m_{j}\omega_{j}=0\right\\}.$ (5.70)

Moreover, $\mathcal{M}$ can be reduced to contain at most $N-1$ integer linear
combinations.

###### Proof.

Combining 5.5 and 5.53 we obtain the formula of the time-averaged
characteristic function of finite-dimensional distribution in the general
case, which is

$\phi_{\theta_{1},\ldots,\theta_{L}}\mathrm{e}^{\mathrm{i}\lambda_{0}X_{0}}\sum_{S\in
G_{N}}\exp\bigg{(}\mathrm{i}\sum_{m_{j}\in
S}m_{j}\tilde{\Theta}_{j}\bigg{)}\prod_{m_{j}\in
S}\mathrm{i}^{m_{j}}J_{m_{j}}(\lambda_{j}R_{j});$ (5.71)

here $G_{N}$ are all $N$-element subsets of integers $m_{j}$ for which
$\sum_{j=1}^{N}m_{j}\omega_{j}=0$. What is left is to show that the above
quantity equals the conditional ensemble-average characteristic function

$\displaystyle\mathbb{E}\left[\mathrm{e}^{\mathrm{i}\left(\sum_{l=1}^{L}\theta_{l}X_{\text{erg}}(t+t_{l})+\sum_{j=1}^{N}\lambda_{j}R_{j}\cos(\tilde{\Theta}_{j}+\omega_{j}t)+\lambda_{0}X_{0}\right)}|X_{0},\\{R_{k}\\},\mathcal{M}\right]$
$\displaystyle=\phi_{\theta_{1},\ldots,\theta_{L}}\mathrm{e}^{\mathrm{i}\lambda_{0}X_{0}}\mathbb{E}\left[\mathrm{e}^{\mathrm{i}\left(\sum_{j=1}^{N}\lambda_{j}R_{j}\cos(\tilde{\Theta}_{j}+\omega_{j}t)\right)}|\\{R_{k}\\},\mathcal{M}\right]$
(5.72)

The factor
$\phi_{\theta_{1},\ldots,\theta_{L}}\mathrm{e}^{\mathrm{i}\lambda_{0}X_{0}}$
already agrees, so we will omit it later on. Next, we expand the remaining
expected value using the Jacobi-Anger identity 5.51, obtaining

$\sum_{S\in M_{N}}\prod_{m_{j}\in
S}\mathrm{i}^{m_{j}}J_{m_{j}}(\lambda_{j}R_{j})\mathbb{E}\left[\mathrm{e}^{\mathrm{i}\sum_{m_{k}\in
S}m_{k}(\tilde{\Theta}_{k}+\omega_{k}t)}|\mathcal{M}\right].$ (5.73)

Because $\theta\mapsto\mathrm{e}^{\mathrm{i}\theta}$ is an injection on
$[0,2\pi)$ and $\tilde{\Theta}_{k}$ differ from $\Theta_{k}$ only by
deterministic constants, the $\sigma$-algebra $\sigma(\mathcal{M})$ is
equivalent to $\sigma(\tilde{\mathcal{M}})$ generated by

$\tilde{\mathcal{M}}=\left\\{\exp\bigg{(}\sum_{m_{j}}m_{j}\tilde{\Theta}_{j}\bigg{)}\colon\sum_{j}m_{j}\omega_{j}=0\right\\}.$
(5.74)

For terms with $S\in G_{N}\subset M_{N}$ the random phases in 5.73 are
$\tilde{\mathcal{M}}$-measurable, therefore also $\mathcal{M}$-measurable,
moreover they agree with the corresponding terms in 5.71. What is left is to
show that the expected value of the remaining terms for $S\notin G_{N}$ is
zero.

Now, for any incommensurable $\omega_{j}$ the corresponding $m_{j}=0$.
Commensurable $\omega_{j}$’s can be divided into subsets of jointly
commensurable numbers, i.e. into blocks $\\{\omega_{k_{i}}\\}$ for which
$\omega_{k_{i}}=\alpha q_{k_{i}}/p_{k_{i}},q_{k_{i}}\in\mathbb{Z},p_{k_{i}}$
and $q_{k_{i}}$ are coprime; different blocks have different incommensurable
factors $\alpha$. Each such block corresponds to a different subset of
independent $\tilde{\Theta}_{j}$, therefore they can be considered separately.

Let us choose one such block and for simplicity of notation, change indices
such that these are $\\{\omega_{1},\omega_{2}\ldots,\omega_{r}\\}$. The
condition $\sum_{j=1}^{r}m_{j}\omega_{j}=0$ is equivalent to the condition
$\sum_{j=1}^{r}m_{j}\eta_{j}=0$, where $\eta_{j}$ are relatively prime
integers obtained by multiplying $\omega_{j}$ by the least common multiple of
$p_{j}$’s. The equation

$\sum_{j=1}^{r}m_{j}\eta_{j}=0$ (5.75)

has exactly $r-1$ linearly independent solutions in integers (Kimme, 1965).
For our one chosen block let us name these solutions
$\\{m_{j}^{1}\\},\\{m_{j}^{2}\\},\ldots,\\{m_{j}^{r-1}\\}$. Any other solution
is a linear combination of the elementary solutions

$m_{j}=\sum_{\rho=1}^{r-1}\nu_{\rho}m_{j}^{\rho},\quad\nu_{j}\in.\mathbb{Z}$
(5.76)

Therefore

$\sum_{j=1}^{r}m_{j}\tilde{\Theta}_{j}=\sum_{j=1}^{r}\sum_{\rho=1}^{r-1}\nu_{\rho}m_{j}^{\rho}\tilde{\Theta}_{j}=\sum_{\rho=1}^{r-1}\nu_{\rho}\sum_{j=1}^{r}m_{j}^{\rho}\tilde{\Theta}_{j},$
(5.77)

and for each block $\mathcal{M}$ depends actually only on $r-1$ variables

$\Xi_{\rho}=\sum_{j=1}^{r}m_{j}^{\rho}\tilde{\Theta}_{j}\ (\mathrm{mod}\
2\pi),$ (5.78)

the rest of the variables are linear combinations of the elementary ones. For
all blocks together $\mathcal{M}$ depends on at most $N-1$ such variables.

The factor in the studied conditional expectation corresponding to the chosen
block is

$\mathbb{E}\left[\exp\bigg{(}\mathrm{i}\sum_{j=1}^{r}m_{j}\tilde{\Theta}_{j}\bigg{)}|\big{\\{}\Xi_{\rho}\big{\\}}_{\rho=1}^{r-1}\right],\quad
m_{j}\in S.$ (5.79)

Because $m_{j}\in S\notin G_{N}$ the sum
$\sum_{j=1}^{r}m_{j}\tilde{\Theta}_{j}\ (\mathrm{mod}\ 2\pi)$ is linearly
independent from the set $\\{\Xi_{\rho}\\}_{\rho=1}^{r-1}$. We prove in the
lemma below that it implies that this sum is also probabilistically
independent from $\\{\Xi_{\rho}\\}_{\rho=1}^{r-1}$ and has uniform
distribution on $[0,2\pi)$. Therefore

$\displaystyle\mathbb{E}\left[\exp\bigg{(}\mathrm{i}\sum_{j=1}^{r}m_{j}\tilde{\Theta}_{j}\bigg{)}|\big{\\{}\Xi_{\rho}\big{\\}}_{\rho=1}^{r-1}\right]=\mathbb{E}\left[\exp\bigg{(}\mathrm{i}\sum_{j=1}^{r}m_{j}\tilde{\Theta}_{j}\bigg{)}\right]$
$\displaystyle=\mathbb{E}\left[\mathrm{e}^{\mathrm{i}\Theta^{\prime}}\right]=\frac{1}{2\pi}\int_{0}^{2\pi}\mathrm{d}\theta\
e^{\mathrm{i}\theta}=0,\quad\Theta^{\prime}\overset{d}{=}\mathcal{U}(0,2\pi).$
(5.80)

We have proven that all elements in the sum 5.73 which contain combinations of
$\Theta_{j}$ linearly independent from elements of $\mathcal{M}$ are zero.
Only $\mathcal{M}$-dependent elements remain, which exactly agrees with the
time-average characteristic function 5.71. This concludes the proof. ∎

######  Lemma 5.5.2.

For i.i.d.
$\\{\Theta_{j}\\}_{j=1}^{N},\Theta_{j}\overset{d}{=}\mathcal{U}(0,2\pi)$ any
$N$ linearly independent integer combinations

$\Xi_{i}=\sum_{j=1}^{N}m_{ij}\Theta_{j}\ (\mathrm{mod}\ 2\pi),\quad
m_{ij}\in\mathbb{Z}$ (5.81)

are a set of jointly independent random variables with distribution
$\Xi_{i}\overset{d}{=}\mathcal{U}(0,2\pi)$.

###### Proof.

Because we work in $\ (\mathrm{mod}\ 2\pi)$ arithmetic, all variables can be
considered to have values in the torus $\mathbb{T}=\mathbb{R}/2\pi\mathbb{Z}$,
$\mathcal{U}(0,2\pi)\equiv\mathcal{U}(\mathbb{T})$. The continuous dual of
$\mathbb{T}^{N}$ is $\mathbb{Z}^{N}$ and this a natural space of parameters of
the characteristic function of the vector $(\Xi_{i})_{i=1}^{N}$.

For a uniform distribution on torus the characteristic function has a very
simple form: it is the Kronecker delta

$\mathbb{E}\left[\mathrm{e}^{\mathrm{i}k\Theta_{j}}\right]=\delta_{k},\quad
k\in\mathbb{Z}.$ (5.82)

This is clear if we think about characteristic function as a Fourier series of
density $1/(2\pi)$ on $\mathbb{T}$. We will show that the multidimensional
characteristic function of $(\Xi_{i})_{i=1}^{N}$ is the product
$\prod_{j=1}^{N}\delta_{k_{j}}$ which corresponds to the distribution
$\mathcal{U}(\mathbb{T}^{N})$ of $(\Xi_{i})_{i=1}^{N}$.

Let us choose any $k_{1},\ldots,k_{N}\in\mathbb{Z}$ and consider
$\sum_{i}k_{i}\Xi_{i}$. We calculate the characteristic function

$\displaystyle\mathbb{E}\left[\exp\bigg{(}\mathrm{i}\sum_{i=1}^{N}k_{i}\Xi_{i}\bigg{)}\right]=\mathbb{E}\left[\exp\bigg{(}\mathrm{i}\sum_{i=1}^{N}k_{i}\sum_{j=1}^{N}m_{ij}\Theta_{j}\bigg{)}\right]$
$\displaystyle=\mathbb{E}\left[\exp\bigg{(}\mathrm{i}\sum_{j=1}^{N}\Theta_{j}\sum_{i=1}^{N}k_{i}m_{ij}\bigg{)}\right]=\prod_{j=1}^{N}\mathbb{E}\left[\exp\Big{(}\mathrm{i}\Theta_{j}\sum_{i=1}^{N}k_{i}m_{ij}\Big{)}\right]$
$\displaystyle=\prod_{j=1}^{N}\delta_{\sum_{i=1}^{N}k_{i}m_{ij}}.$ (5.83)

The above product equals 1 and not 0 if, and only if for all $j$ we have
$\sum_{i=1}^{N}k_{i}m_{ij}=0$. But, the linear integer combinations 5.81 are
linearly independent which is equivalent to saying that this is true if, and
only if $k_{1}=k_{2}=\ldots=k_{N}=0$. The last product is
$\prod_{j=1}^{N}\delta_{k_{j}}$ which was to be proven.

∎

This general theorem completely determines the behaviour of time-averaged
observables for stationary Gaussian processes and has important practical
applications, allowing for the statistical analysis of non-ergodic stationary
Gaussian models.

As an example, let us come back to the time-averaged characteristic function
5.69. For the case of a process with two spectral points with frequencies
$\omega_{1}=1,\omega_{2}=2$, the numerically calculated time-averages are
shown shown as yellow lines in Fig. 5.4. The integer combinations in the
family $\mathcal{M}$ from Eq. 5.70 are exactly $m_{1}=2,m_{2}=-1$ and
multiples of these. Therefore the time-averaged characteristic function
depends only on $2\Theta_{1}-\Theta_{2}\ (\mathrm{mod}\ 2\pi)$. Indeed, any
yellow line in Fig. 5.4 corresponds to many different random choices of
$\Theta_{1},\Theta_{2}$, and these lines can be parametrised only by the
number $c\in[0,2\pi)$ defined as $2\Theta_{1}-\Theta_{2}=c\ (\mathrm{mod}\
2\pi)$. Similarly, green lines on the same figure depend only on
$3\Theta_{1}-2\Theta_{2}\ (\mathrm{mod}\ 2\pi)$, because we have
$m_{1}=3,m_{2}=-2$, and so on.

It may seem counter-intuitive that the rationality or irrationality of the
number $\omega_{2}$, which cannot be experimentally studied, affects the
numerical simulations and the behaviour of the real systems. This apparent
paradox disappears if we carefully analyse the behaviour of the time-average
in the two essential cases

  * •

For $\omega_{2}=p/q$ with coprime $p,q$ the trajectory will have period $2\pi
q$. For $q$ sufficiently larger than the experimental time ($q\gg T$), the
periodicity will be unobservable during the measurement and the time-averaged
observables will seem independent from initial phases.

  * •

For irrational $\omega_{2}=p/q+\epsilon$ with coprime $p,q$ and $\epsilon q\ll
1$ the time-averaged observables will not depend on initial phases in the
long-time limit $T\to\infty$, however the process will be very close to a
periodic one, therefore the convergence will be slow.

So, we realise that in the real case with finite time of the experiment
$T<\infty$ the practically significant property is how close $\omega_{2}$ is
to a irreducible fraction with a small denominator. Analogously, for multiple
$\omega_{k}$ we only need to determine how close they are to a set of
commensurable numbers with simple rational ratios.

In order to illustrate this fact numerically we simulated the process
$R_{1}\cos(\Theta_{1}+t)+R_{2}\cos(\Theta_{2}+\omega_{2}t)$ with fixed
$R_{1}=2,R_{2}=3$ and calculated the variance of the time-averaged
characteristic function at one point $\overline{\phi}(2.5)$. For irrational
$\omega_{2}$ and in the limit $T\to\infty$ this quantity should be equal to
zero, however in the numerical experiment it is a positive function of
$\omega_{2}$ which indeed measures how close $\omega_{2}$ is to a simple
irreducible fraction. See Fig. 5.5 where the peaks of the variance indicate
the positions of the simplest irreducible fractions like $1/2,1/3,2/3,1/4,3/4$
_etc_.

Figure 5.5: Numerical estimation of the variance
$\operatorname*{{var}}[\overline{\phi}(2.5)]$ obtained using $10^{3}$ samples
of trajectories with lengtht $T=200$ and values of $\omega_{2}$ taken as one
thousand uniformly scattered numbers in the interval $[0,1]$ stored in format
double.

Another result, which may be somehow unexpected, is that according to the
calculation in Section 5.3 even for commensurable $\\{\omega_{k}\\}_{k}$ the
time-average second-order properties depend only on $\\{R_{k}\\}_{k}$ and the
dependence of initial phases $\\{\Theta_{k}\\}_{k}$ is lost.

######  Proposition 5.5.1.

For any stationary Gaussian process $X$, the time-average covariance structure
is the ensemble-average structure conditioned by the $\sigma$-algebra
$\sigma(X_{0},\\{R_{k}\\}_{k})$, where $R_{k}$ are the amplitudes of the atoms
of the spectral measure and $X_{0}$ is the constant term. The result is true
even for rationally commensurable frequencies $\\{\omega_{k}\\}$ of spectral
points and infinite number of them, $N=\infty$.

###### Proof.

We begin with calculating the conditional ensemble-average covariance. The
conditional mean equals

$\mathbb{E}[X(t)|X_{0},\\{R_{k}\\}]=X_{0},\quad\text{as
}\mathbb{E}[\cos(\Theta_{k})]=0\text{, and }\mathbb{E}[X_{\text{erg}}(t)]=0$
(5.84)

Next we fix $t_{1},t_{2}$ and use the independence of the $X_{\text{erg}}$ and
$X_{k}$’s.

$\displaystyle\
\mathbb{E}\big{[}(X(t_{1})-X_{0})(X(t_{2})-X_{0})|X_{0},\\{R_{k}\\}\big{]}$
$\displaystyle=\mathbb{E}\big{[}X_{\text{erg}}(t_{1})X_{\text{erg}}(t_{2})|X_{0},\\{R_{k}\\}\big{]}+\sum_{k=1}^{\infty}\mathbb{E}\big{[}X_{k}(t_{1})X_{k}(t_{2})|X_{0},\\{R_{k}\\}\big{]}$
$\displaystyle=\mathbb{E}\big{[}X_{\text{erg}}(t_{1})X_{\text{erg}}(t_{2})\big{]}+\sum_{k=1}^{\infty}R_{k}^{2}\mathbb{E}\big{[}\cos(\Theta_{k}+\omega_{k}t_{1})\cos(\Theta_{k}+\omega_{k}t_{2})\big{]}$
$\displaystyle=r(t_{2}-t_{1})+\sum_{k=1}^{\infty}R_{k}^{2}\frac{1}{2\pi}\int_{0}^{2\pi}\mathrm{d}x\
\cos(x+\omega_{k}t_{1})\cos(x+\omega_{k}t_{2})$
$\displaystyle=r(t_{2}-t_{1})+\frac{1}{2}\sum_{k=1}^{\infty}R_{k}^{2}\cos(\omega_{k}(t_{2}-t_{1})).$
(5.85)

The time-average consist of an integral from the parts
$X_{\text{erg}}(t_{i}+\tau)X_{\text{erg}}(t_{j}+\tau)$,$X_{k}(t_{i}+\tau)X_{k}(t_{j}+\tau)$,
$X_{k_{1}}(t_{i}+\tau)X_{k_{2}}(t_{j}+\tau),k_{1}\neq k_{2}$ and
$X_{\text{erg}}(t_{i}+\tau)X_{k}(t_{j}+\tau)$, $i,j\in\\{1,2\\}$; we call them
$I_{1},I_{2},I_{3},I_{4}$, respectively. All sums are absolutely convergent:
we can commute summation, integration and taking limit $T\to\infty$.

Time-average $T^{-1}I_{1}$ converges to $r(t_{2}-t_{1})$ because
$X_{\text{erg}}$ is ergodic. For $T^{-1}I_{2}$ we have

$\displaystyle\frac{1}{T}I_{4}=\frac{1}{T}R_{k}^{2}\int_{0}^{T}\mathrm{d}\tau\
\cos(\tau+\omega_{k}t_{i})\cos(\tau+\Theta_{k}+\omega_{k}t_{j})$
$\displaystyle=\frac{1}{2}R_{k}^{2}\cos(\omega_{k}(t_{2}-t_{1}))$
$\displaystyle+\frac{1}{\omega_{k}T}\cos(2\Theta_{k}+\omega_{k}(t_{1}+t_{2}+T))\sin(\omega_{k}T)\xrightarrow{T\to\infty}\frac{1}{2}R_{k}^{2}\cos(\omega_{k}(t_{2}-t_{1})).$
(5.86)

Therefore we need to prove that $T^{-1}I_{3}$ and $T^{-1}I_{4}$ decay to 0.
For $T^{-1}I_{3}$ it is straightforward, denoting
$\omega_{\pm}\mathrel{\mathop{:}}=\omega_{k_{2}}\pm\omega_{k_{1}},\Theta_{\pm}\mathrel{\mathop{:}}=\Theta_{k_{2}}\pm\Theta_{k_{1}}$
we get

$\displaystyle\
\frac{1}{T}I_{3}=\frac{1}{2T\omega_{-}}\big{(}\sin(\omega_{k_{1}}t_{i}-\omega_{k_{2}}t_{j}-\Theta_{-})-\sin(\omega_{k_{1}}t_{i}-\omega_{k_{2}}t_{j}-\Theta_{-}-T\omega_{-})\big{)}$
$\displaystyle+\frac{1}{2T\omega_{+}}\big{(}-\sin(\omega_{k_{1}}t_{i}+\omega_{k_{2}}t_{j}+\Theta_{+})-\sin(\omega_{k_{1}}t_{i}+\omega_{k_{2}}t_{j}+\Theta_{+}+T\omega_{+})\big{)}\xrightarrow{T\to\infty}0.$
(5.87)

As for $T^{-1}I_{4}$, we will show that the time-average of
$X_{\text{erg}}(t_{i}+\tau)R_{k}\exp(\mathrm{i}\Theta_{k}+\mathrm{i}\omega_{k}\tau)$
converges, which is an equivalent condition. The factor
$R_{k}\exp(\mathrm{i}\Theta_{k})$ can be brought outside the integral,
therefore only showing the convergence of time-average of
$X_{\text{erg}}(t_{i}+\tau)\exp(\mathrm{i}\omega_{k}\tau)$ is required. The
latter is

$\frac{1}{T}\int_{0}^{T}\mathrm{d}\tau
X_{\text{erg}}(\tau)\mathrm{e}^{\mathrm{i}\omega_{k}\tau}.$ (5.88)

The limit $T\to\infty$ of the above formula exists almost surely, the argument
is the same as at the beginning of the lemma. We will prove it is 0. Let us
calculate $\mathbb{E}|\boldsymbol{\cdot}|^{2}$ of 5.88

$\displaystyle\
\frac{1}{T^{2}}\int_{0}^{T}\mathrm{d}\tau_{1}\int_{0}^{T}\mathrm{d}\tau_{2}\
r(\tau_{2}-\tau_{1})\mathrm{e}^{\mathrm{i}\omega_{k}(\tau_{2}-\tau_{1})}$
$\displaystyle=\frac{1}{T^{2}}\int_{0}^{T}\mathrm{d}\tau_{1}\int_{0}^{T}\mathrm{d}\tau_{2}\int_{\mathbb{R}}\mathrm{d}\sigma(\omega)\mathrm{e}^{\mathrm{i}\omega(\tau_{2}-\tau_{1})}\mathrm{e}^{\mathrm{i}\omega_{k}(\tau_{2}-\tau_{1})}$
$\displaystyle=\frac{1}{T^{2}}\int_{\mathbb{R}}\mathrm{d}\sigma(\omega)\int_{0}^{T}\mathrm{d}\tau_{1}\int_{0}^{T}\mathrm{d}\tau_{2}\mathrm{e}^{\mathrm{i}(\omega+\omega_{k})(\tau_{2}-\tau_{1})}$
$\displaystyle=2\int_{\mathbb{R}}\mathrm{d}\sigma(\omega-\omega_{k})\frac{1-\cos(\omega
T)}{(\omega
T)^{2}}\xrightarrow{\mathrm{T}\to\infty}\sigma(\\{\omega_{k}\\})=0.$

That shows the mean-square and almost sure convergence to 0. ∎

It means that the additional memory structure induced by the commensurability
of $\\{\omega_{k}\\}_{k}$ is purely non-linear and it is possible to detect it
only using higher-order statistics, e.g. the time-averaged characteristic
function. This fact do not contradict the purely linear dependence structure
of the Gaussian processes, because it applies only to the time-averages, which
in the non-ergodic case make use only of a part of the full information
contained in the process.

The exemplary covariance estimation is shown in Fig. 5.6, where we sampled the
stationary Ornstein-Uhlenbeck process with mean-returning parameter
$\lambda=3$, and the addition of two spectral components $2\cos(t)$ and
$\cos(2t)$. The covariance conditioned by $R_{1}=2,R_{2}=1$ is

$r(t)=\mathrm{e}^{-\lambda t}+2\cos(t)+\frac{1}{2}\cos(2t)$ (5.89)

and it agrees with the estimated time-averaged covariance shown in Fig. 5.6.

Figure 5.6: Comparison between estimated time-average covariance function and
a theoretical ensemble one for a stationary non-ergodic process with
commensurable frequencies of the spectral points.

###  5.6 Ergodicity of Langevin equations

In this section we will make use of three basic facts.

######  Proposition 5.6.1.

Let us consider any finite measure $\sigma$ and a measurable function $f$,
defined $\sigma$-almost everywhere. In this case:

  * i)

If $\sigma$ is absolutely continuous, then $\sigma f$ is absolutely
continuous.

  * ii)

If $\sigma$ is Rajchmann, then $\sigma f$ is Rajchmann; equivalently: if
$\widehat{\sigma}$ decays, then $\widehat{\sigma f}$ also decays.

  * iii)

If $\sigma$ is continuous, then $\sigma f$ is continuous.

The function $f$ is often called the spectral gain.

###### Proof.

Fact i) fallows directly from the definition of $\sigma f$. The proof of Fact
iii) is also simple: $\sigma(\\{x_{0}\\})f(x_{0})\neq 0$ only if
$\sigma(\\{x_{0}\\})\neq 0$.

For Fact ii) the argument is a little more complicated. Because the measure
$\sigma$ is finite, trigonometric functions are dense in
$L^{1}(\mathbb{R},\sigma)$. We take a sequence

$f_{k}(\omega)=\sum_{j=1}^{N_{k}}A_{j}\mathrm{e}^{\mathrm{i}t_{j}\omega},$
(5.90)

such that $\|f-f_{k}\|_{L^{1}(\mathbb{R},\sigma)}\to 0$. Note that

$\widehat{\sigma
f_{k}}(t)=\sum_{j=1}^{N_{k}}A_{j}\widehat{\sigma}(\omega+\omega_{j})\xrightarrow{t\to\infty}0,$
(5.91)

in other words $\sigma f_{k}$ are Rajchamann. Additionally

$\big{|}\widehat{\sigma f}(t)-\widehat{\sigma
f_{k}}(t)\big{|}=\left|\int_{\mathbb{R}}\sigma(\mathrm{d}\omega)\
\mathrm{e}^{\mathrm{i}\omega
t}\big{(}f(\omega)-f_{k}(\omega)\big{)}\right|\leq\|f-f_{k}\|_{L^{1}(\mathbb{R},\sigma)},$
(5.92)

so sequence of functions $\widehat{\omega f_{k}}$ converge uniformly to
$\widehat{\omega f}$. Therefore we can commute the limits
$k\to\infty,t\to\infty$ and obtain

$\lim_{t\to\infty}\widehat{\sigma f}(t)=0.$ (5.93)

∎

Considered together, these three facts guarantee that if a Gaussian process
$X$ has spectral measure $\sigma$, then the Gaussian process $Y$ with spectral
measure $\sigma f$ inherits all ergodic properties (ergodicity, mixing) from
$X$. In fact, in some cases the time-averages of $Y$ can behave more regularly
than those of $X$, which will be shown later in this section.

The above proposition can be used to determine the ergodic behaviour of
various transformations of a given process. The process of type

$Y(t)=\sum_{k}a_{k}X(t-T_{k}),$ (5.94)

discussed in Chapter 2, appears directly e.g. in the biological systems in
which delayed responses are to be expected (Barbuti et al., 2011).

Proposition 2.2.3 shows that $Y$ has a spectral measure of $X$ multiplied by a
continuous function, so it inherits the ergodic properties of the process $X$.
The case

$Y(t)=\sum_{k}a_{k}\frac{\mathrm{d}^{k}}{\mathrm{d}t^{k}}X(t),$ (5.95)

considered in Proposition 2.2.4 is similar. Such a $Y$ inherits the ergodic
properties of $X$, but can also be ergodic whereas $X$ is not. The function
$\omega\mapsto\left|\sum_{k}a_{k}(\mathrm{i}\omega)^{k}\right|^{2}$ may have
zeros and if such zero agrees with the position of the spectral point of $X$,
the process $Y$ does not contain this spectral point. The simplest such case
is when $X$ has exactly one spectral point at $\omega=0$, that is it contains
a time-independent Gaussian constant $X_{0}$. In this situation any time-
averaged observable which depends on the mean of $X$ does not converge to the
ensemble-average. However, the derivative
$\frac{\mathrm{d}}{\mathrm{d}t}X(t)$, corresponding to the spectral gain
function $\omega\mapsto\omega^{2}$, does not contain $X_{0}$ and is ergodic.

For $Y(t)=g*X(t)$ from Proposition 2.2.5 the result is also analogous, $Y$
inherits the ergodic properties of $X$. If the zeros of $\widehat{g}$ agree
with the spectral points of $X$, the process $Y$ may be ergodic whereas $X$ is
not. Moreover, if the singular non-mixing measure of $X$ is contained in the
domain outside of $\widehat{Y}$ support, $Y$ may be mixing when $X$ is non-
mixing.

Special care must be taken in the case when we consider the convolution with
$g$ which is not rapidly decaying, e.g. $g\in L^{2}(\mathbb{R})$ but $g\notin
L^{1}(\mathbb{R})$. The Fourier transform of such $g$ has jumps. The well-
known result from Fourier theory (Champeney, 1987) states that if
$\widehat{g}$ has a jump at $\omega_{0}$, then $\int_{-T}^{T}\mathrm{d}t\
g(t)\mathrm{e}^{\mathrm{i}\omega_{0}t}$ converges as $T\to\infty$ to the value
$(\widehat{g}(\omega_{0}^{-})+\widehat{g}(\omega_{0}^{+}))/2$, i.e. to the
exact middle of the discontinuity of $\widehat{g}$. If the process $X$ has a
spectral atom at the exact frequency $\omega_{0}$, then for the process $Y$
this spectral point will be modulated by
$|\widehat{g}(\omega_{0}^{-})+\widehat{g}(\omega_{0}^{+})|^{2}/4$ as long as
the filter is applied symmetrically as a limit of convolutions with functions
$g$ which are supported on interval $[-T,T]$.

The harmonic representation of the convolution vastly increases the number of
models for which it is easy to study ergodicity using Fourier methods, as
using convolution is one of the most often chosen methods to model time-
invariant linear responses of the system (see also next Sec. 5.6). One of the
most common examples that appear in practice is $g$ being a one or two sided
exponent decaying with rate $\lambda$. These choices correspond to the
spectral responses $1/(\lambda^{2}+\omega^{2})$ or
$4\lambda^{2}/(\lambda^{2}+\omega^{2})$, respectively.

One other practical consequence is that one can filter out non-ergodicity from
the data. Using estimators of power spectral density (e.g. periodogram
(Brockwell and Davis, 2006)) the locations of spectral points can be estimated
and the corresponding non-ergodicity removed by using any filter with zeros
its spectral gain function at their frequencies. The simplest choice of such
filter is the smoothing

$\tilde{X}(t)=\int\limits_{t-\pi/\omega_{0}}^{t+\pi/\omega_{0}}\mathrm{d}s\
X(s),$ (5.96)

which integrates the spectral component $\omega_{0}$ over its period,
therefore removing it. It corresponds to the filter $g_{\omega_{0}}$ and
spectral gain $\widehat{g}_{\omega_{0}}$

$g_{\omega_{0}}(t)=\begin{cases}1,&|t|\leq\pi/\omega_{0},\\\
0,&|t|>\pi/\omega_{0},\end{cases}\quad\widehat{g}_{\omega_{0}}(\omega)=\sqrt{\frac{2\omega_{0}^{2}}{\pi^{3}}}\operatorname{sinc}\left(\frac{\omega\pi}{\omega_{0}}\right).$
(5.97)

For a multiple number of spectral points one may use the filter
$g_{\omega_{1}}*g_{\omega_{2}}\ldots g_{\omega_{N}}$ or any other with a
suitable spectral gain. The gain $|\widehat{g}|^{2}$ behaves like
$(\omega-\omega_{0})^{2}$ near $\omega_{0}$, therefore it removes the spectral
point $\omega_{0}$ in a numerically stable manner. If one is more interested
in sure removal of the non-ergodicity near the location $\omega_{0}$ than not
distorting the spectrum, one can use a spectral gain which is more flat around
$\omega_{0}$, e.g. using a triangular function filter guarantees the
asymptotic behaviour $(\omega-\omega_{0})^{4}$. An other useful choice is
spectral gain

$\widehat{g}(\omega)=\begin{cases}1-\mathrm{e}^{-(\omega-\omega_{0})^{2}/c},&|\omega|\leq
L;\\\ 0,&|\omega|>L,\end{cases}$ (5.98)

which allows for the calibration of the level of distortion around
$\omega_{0}$ (parameter $c$) and frequency cut-off (parameter $L$). The
corresponding filter can be expressed using error, Gaussian and trigonometric
functions, so it can be easily computed for the purpose of statistical usage.

The above approach can be understood as, instead of using original observables
$f(X)$, to use the modified observable $\widetilde{f}(X)=f(Y)$ for which the
time- and ensemble-averages coincide, even when this is not generally the
case. Analysing the transformed process $Y$ instead of $X$ may be more
difficult, as the properties of $X$ are distorted by filtering, albeit in a
controlled manner. However, for a small number of spectral points it is
manageable, moreover it can be used as an effective method of localising
spectral points: if the observables of the filtered process behave like it is
ergodic, it is a statistical verification of a good choice of these locations.
Further analysis can be performed on the filtered process $Y$, which is
ergodic, or by staying with the original $X$ and using methods from Section
5.5.

Identical reasoning applies to stochastic differential equations. The simplest
case of the Langevin equation

$\mathrm{d}X=-\mathrm{d}t\lambda X+\mathrm{d}B$ (5.99)

is solved by the Ornstein-Uhlenbeck process, as discussed in Section 3.4. It
has a spectral density $s_{X}(\omega)=(\lambda^{2}+\omega^{2})^{-1}$ so it is
clearly mixing and that would be the case also for equation with more general
mixing force than $\mathrm{d}B$, e.g. $\mathrm{d}B_{H}$. This result is much
more general. Any stationary solution of the linear differential system

$a_{k}\frac{\mathrm{d}^{k}}{\mathrm{d}t^{k}}X(t)+a_{k-1}\frac{\mathrm{d}^{k-1}}{\mathrm{d}t^{k-1}}X(t)+\ldots+a_{0}X(t)=F(t),$
(5.100)

has a casual solution given by a proper Green’s function and spectral
representation

$X(t)=\int_{\mathbb{R}}\mathrm{d}S_{F}(\omega)\
\frac{1}{\sum_{j=1}^{k}a_{j}(\mathrm{i}\omega)^{j}}\mathrm{e}^{\mathrm{i}\omega
t},$ (5.101)

It is clear that $X$ the inherits ergodic properties of $F$. In contrast to
the case when $X$ was an effect of a linear filter on $Y$, in this case $X$
cannot even be ergodic if $F$ is not, as the rational function
$1/|\sum_{j}a_{j}(\mathrm{i}\omega)^{j}|^{2}$ does not have zeros.

Very similar reasoning applies to systems with richer memory structure, the
most prominent being our main subject of interest, the generalized Langevin
equation, which we will here analyse in its full form

$m\frac{\mathrm{d}^{2}}{\mathrm{d}t^{2}}X(t)=-\kappa
X-\int_{-\infty}^{t}\\!\\!\\!\mathrm{d}\tau\
\frac{\mathrm{d}}{\mathrm{d}t}X(\tau)K(t-\tau)+F(t).$ (5.102)

Its stationary solution has the form

$X(t)=\int_{\mathbb{R}}\mathrm{d}S_{F}(\omega)\
\frac{1}{-m\omega^{2}+\kappa+\mathrm{i}\omega\widehat{K}(\omega)},$ (5.103)

where $\mathrm{d}S_{F}$ is the spectral process of the force $F$. The GLE with
this form has a visible similarity to the Newton equation and can model
subdiffusion in the case when $F$ is fractional Gaussian noise (Viñales and
Despósito, 2006; Kou, 2008), and the kernel is $K(t)=t^{-\alpha},0<\alpha<2$.
Such a kernel is outside $L^{2}$, but the harmonic representation above is
still valid if we interpret $\widehat{K}$ in the generalized sense, as
commented in Section 2.3. Using this approach we can use tables of generalized
Fourier transform to obtain
$\widehat{K}(\omega)=|\omega|^{\alpha-1}\Gamma(1-\alpha)\mathrm{e}^{\mathrm{i}\,\text{sgn}(\omega)\pi(\alpha-1)/2}$
(Champeney, 1987).

For $\kappa>0$ the fraction in Eq. 5.103 behaves like a constant near $0$ and
like $\omega^{-2}$ in the limit $\omega\to\pm\infty$. Therefore the spectral
measure $\sigma_{X}$ is generally finite. When there is no external potential,
$\kappa=0$, the particle is not confined, so physically speaking there should
be no stationary solution. It is also visible from Eq. 5.103 in the fact that
if there was a stationary solution, its spectral measure $\sigma_{X}$ would
exhibit singularity $\omega^{-2}$ at $0$, which is a contradiction of
$\sigma_{X}$ being a finite or truncated measure. In this case however the
velocity would be a stationary process and $X$ could be represented as

$X(t)=\int_{\mathbb{R}}\mathrm{d}S_{F}\frac{\mathrm{e}^{\mathrm{i}\omega
t}-1}{-m\omega^{2}+\mathrm{i}\omega\widehat{K}(\omega))}.$ (5.104)

In both cases the position $X$ or the velocity
$V=\frac{\mathrm{d}}{\mathrm{d}t}X$ inherit the ergodic properties of $F$.

The Fourier space approach gives also additional insight into the physical
origin of the non-ergodicity. In Chapter 3 we presented a derivation of the
GLE based on discrete and continuous field baths. For a pure field bath the
spectral measure of $F$ was determined by the Fourier transform of the
coupling function

$\sigma_{F}(\mathrm{d}\omega)=\mathrm{d}\omega\ |\widehat{\rho}|^{2},$ (5.105)

which lead to an immediate conclusion.

######  Proposition 5.6.2.

The stationary solution of the Hamiltonian system from Section 3.1, in which
only the continuous heat bath is present, is mixing.

Conversely, the force related to the discrete heat bath has a pure atomic
spectrum supported on self-frequencies of the bath harmonic oscillators

$\sigma_{F}(\mathrm{d}\omega)=\sum_{k}m_{k}\frac{\gamma_{k}^{2}}{2\omega_{k}^{2}}\big{(}\delta(\mathrm{d}\omega-\omega_{k})+\delta(\mathrm{d}\omega+\omega_{k})\big{)}.$
(5.106)

This fact immediately prohibits the ergodicity of the solution. Actually, if
the fluctuation-dissipation theorem holds, even the stationary solution does
not exist, as the integral $\int_{-\infty}^{t}\\!\\!\\!\mathrm{d}\tau\
X(\tau)K(t-\tau)$ does not converge for stationary $X$ and non-decaying $K$.
The non-stationary solutions of the equation with the finite-time integral
$\int_{0}^{t}$ are still well-defined.

From the point of view of spectral theory, for non-decaying $K$ its Fourier
transform can be defined only as a distribution and the existence of the
stationary solution would require dividing by a distribution, which cannot be