Random coefficient autoregressive processes describe Brownian yet non-Gaussian diffusion in heterogeneous systems

TL;DR

This paper connects physical models of heterogeneous diffusion with statistical autoregressive models, enabling better analysis of complex diffusion behaviors in biological systems.

Contribution

It establishes a theoretical link between Langevin equations with variable coefficients and random coefficient autoregressive processes, bridging physics and statistics.

Findings

Provides a framework to identify Brownian yet non-Gaussian diffusion in data.

Enables application of statistical methods to physical stochastic processes.

Bridges the gap between physical heterogeneity models and time series analysis.

Abstract

Many studies on biological and soft matter systems report the joint presence of a linear mean-squared displacement and a non-Gaussian probability density exhibiting, for instance, exponential or stretched-Gaussian tails. This phenomenon is ascribed to the heterogeneity of the medium and is captured by random parameter models such as "superstatistics" or "diffusing diffusivity". Independently, scientists working in the area of time series analysis and statistics have studied a class of discrete-time processes with similar properties, namely, random coefficient autoregressive models. In this work we try to reconcile these two approaches and thus provide a bridge between physical stochastic processes and autoregressive models. We start from the basic Langevin equation of motion with time-varying damping or diffusion coefficients and establish the link to random coefficient autoregressive processes. By exploring that link we gain access to efficient statistical methods which can help to identify data exhibiting Brownian yet non-Gaussian diffusion.

Equations

$${\delta }_X^2(t):=\mathbb{E}[X(t{)}^2]=2Dt$$

$$p_{X}(x;t)=\frac{1}{\sqrt{4\pi Dt}}\exp\mathopen{}\mathclose{{}\left(-\frac{x^{2}}{4Dt}}\right).$$

$${\delta }_X^2(t)\sim D_{\alpha }{t}^{\alpha }$$

$$p_{X}(x;t)=\frac{1}{\sqrt{2Dt}}\exp\mathopen{}\mathclose{{}\left(-\sqrt{\frac{2}{Dt}}|x|}\right),$$

$$p_{X}(x;t)\propto\exp\mathopen{}\mathclose{{}\left(-\mathopen{}\mathclose{{}\left(\frac{|x|}{ct^{\alpha/2}}}\right)^{\delta}}\right),$$

$$\mathrm{d}X(t)=\sqrt{D}\mathrm{d}B(t),$$

$$\mathrm{d}X(t)=\sqrt{D(t)}\mathrm{d}B(t)\text{or equivalently}X(t)={\int }_0^t\sqrt{D(s)}\mathrm{d}B(s).$$

$$V_k-{\varphi }^1{V}_{k-1}+{\varphi }^2{V}_{k-2}+\dots {\varphi }^p{V}_{k-p}=Z_k+{\theta }^1{Z}_{k-1}+{\theta }^2{Z}_{k-2}+\dots +{\theta }^q{Z}_{k-q}.$$

$$Z_k^d=\sum _{j=1}^{\infty }\frac{\Gamma (j+d)}{\Gamma (d)\Gamma (j+1)}{Z}_{k-j}.$$

$$V_k={\Sigma }_k{Z}_k,{\Sigma }_k^2={\alpha }_0+{\alpha }_1{V}_{k-1}^2+{\alpha }_2{V}_{k-2}^2+\dots +{\alpha }_q{V}_{k-q}^2.$$

$$V_k={\Sigma }_k{Z}_k,{\Sigma }_k^2={\alpha }_0+\sum _{i=1}^q{\alpha }_i{V}_{k-i}^2+\sum _{j=1}^p{\beta }_j{\Sigma }_{k-j}^2.$$

$$V_k-{\Phi }_k^1{V}_{k-1}+{\Phi }_k^2{V}_{k-2}+\dots {\Phi }_k^p{V}_{k-p}={\Theta }_k^0{Z}_k+{\Theta }_k^1{Z}_{k-1}+{\Theta }_k^2{Z}_{k-2}+\dots +{\Theta }_k^q{Z}_{k-q}.$$

$$m\frac{\mathrm{d}V(t)}{\mathrm{d}t}=-\beta V(t)+F(t),$$

$$\mathrm{d}V(t)=-\Lambda (t)V(t)\mathrm{d}t+\sqrt{D(t)}\mathrm{d}B(t).$$

$$V_k-(1-{\Lambda }_k\Delta t)V_{k-1}=\sqrt{D_k}\Delta B_k.$$

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

$$X_k-{\varphi }^1{X}_{k-1}-{\varphi }^2{X}_{k-2}=Z_k+\theta Z_{k-1}$$

$${\varphi }^1$$

$${\varphi }^2$$

$${\text{rme}}^{{\int }_0^{k\Delta t}\Lambda (s)\mathrm{d}s}{V}_k-{\text{rme}}^{{\int }_0^{(k-1)\Delta t}\Lambda (s)\mathrm{d}s}{V}_{k-1}={\int }_{(k-1)\Delta t}^{k\Delta t}{\text{rme}}^{{\int }_0^s\Lambda (s^{\prime })\mathrm{d}{s}^{\prime }}\sqrt{D(s)}\mathrm{d}B(s).$$

$$V_k-{\Phi }_k{V}_{k-1}=Z_k$$

$${\Phi }_k:={\text{rme}}^{-{\int }_{(k-1)\Delta t}^{k\Delta t}\Lambda (s)\mathrm{d}s},Z_k:={\int }_{(k-1)\Delta t}^{k\Delta t}{\text{rme}}^{-{\int }_s^{k\Delta t}\Lambda (s^{\prime })\mathrm{d}{s}^{\prime }}\sqrt{D(s)}\mathrm{d}B(s).$$

$$\mathbb{E}[Z_k^2|\Lambda ,D]$$

$$={\int }_{(k-1)\Delta t}^{k\Delta t}{\text{rme}}^{-2{\int }_s^{k\Delta t}\Lambda (s^{\prime })\mathrm{d}{s}^{\prime }}D(s)\mathrm{d}s.$$

$$V_k-{\Phi }_k{V}_{k-1}={\Theta }_k{W}_k,$$

$${\Phi }_k\approx {\text{rme}}^{-\Delta t{\Lambda }_k},{\Theta }_k\approx \sqrt{\frac{D_k}{2{\Lambda }_k}}\sqrt{1-{\text{rme}}^{-2\Delta t{\Lambda }_k}}\approx \sqrt{\frac{\Delta t{D}_k}{-2\ln ({\Phi }_k)}}\sqrt{1-{\Phi }_k^2}.$$

$$V_k$$

$$+{\Phi }_k\cdots {\Phi }_{k-n+2}{\Theta }_{k-n+1}{W}_{k-n+1}.$$

$$V_{k}=\sum_{i=0}^{\infty}\Big{(}\prod_{j=0}^{i-1}\Phi_{k-j}\Big{)}\Theta_{k-i}W_{k-i}.$$

$$S^{2}\mathrel{\mathop{:}}=\mathbb{E}\mathopen{}\mathclose{{}\left[V_{k}^{2}|\Phi,\Theta}\right]=\sum_{i=0}^{\infty}\Big{(}\prod_{j=0}^{i-1}\Phi_{k-j}^{2}\Big{)}\Theta_{k-i}^{2}.$$

$$\kappa_{V}\mathrel{\mathop{:}}=\frac{\mathbb{E}\mathopen{}\mathclose{{}\left[V_{k}^{4}}\right]}{\mathbb{E}\mathopen{}\mathclose{{}\left[V_{k}^{2}}\right]^{2}}-3=3\frac{\mathbb{E}\mathopen{}\mathclose{{}\left[(S^{2}-\mathbb{E}\mathopen{}\mathclose{{}\left[S^{2}}\right])^{2}}\right]}{\mathbb{E}\mathopen{}\mathclose{{}\left[S^{2}}\right]^{2}}$$

$$S_{\alpha }^2=\frac{{\Theta }_{\alpha }^2}{1-{\Phi }_{\alpha }^2},$$

Full text

1 Introduction

Brownian motion, one of the most fundamental processes in non-equilibrium statistical physics, describes the motion of a passive colloidal particle in a thermal fluid environment. It was observed even in ancient times, for instance, by Roman philosopher Lucretius [1]. The modern scientific interest was started by botanist Robert Brown [2]; it was later studied theoretically by Einstein, Sutherland, Smoluchowski, and Langevin between 1905 and 1908 [3, 4, 5, 6]. Two fundamental properties are typically associated with Brownian motion, namely, the linear growth

[TABLE]

in time of the mean-squared displacement (MSD) with the diffusion coefficient $D$, and the Gaussian probability density function (PDF)

[TABLE]

Alternatively to the MSD (1) Brownian motion is characterised by a $1/f^{2}$ frequency dependence of the associated ensemble and single trajectory power spectra [7, 8].

Deviations from the linear time dependence (1) of the MSD are observed routinely in a large variety of systems, in particular, in the power-law form

[TABLE]

with the anomalous diffusion coefficient $D_{\alpha}$ [9, 10, 11]. In biological cells or other complex liquids both subdiffusion with $0<\alpha<1$ [12, 13, 14, 15, 16, 17] and superdiffusion with $1<\alpha<2$ [18, 19, 20] are measured, see the recent reviews [21, 22]. Anomalous diffusion are effected when the increments of the stochastic process are no longer independent, when the variance of the step length or the mean-step time diverges, as well as due to the tortuosity of the embedding space. The associated PDF of these processes may have both Gaussian and non-Gaussian shapes [9, 10, 11].

Of late, a new class of diffusive dynamics has come into focus, following numerous reports in soft matter, biological and other complex systems: in these systems the MSD is normal of the form (1) with invariable coefficient $D$, however, the PDF is non-Gaussian and often found to be of the distinct exponential shape ("Laplace distribution")

[TABLE]

see [23, 24, 25, 26, 27] as well as the extensive list of references in [28]. These "Brownian yet non-Gaussian" processes along with a more general class of non-Gaussian PDFs, discussed in more detail below, are in the focus of this study.

Our goal here is, however, different from the previous modelling approaches. Namely, here we try to establish a direct connection of the physical models for Brownian yet non-Gaussian diffusion and a class of processes ubiquitously used in time series analysis, the so-called autoregressive models with random coefficients [29, 30, 31].

Thus, we introduce the time series methods to the modeling of non-Gaussian diffusion, which enables us to provide information about the distribution of the process, showing the influence of some effective properties of the heterogeneous medium. In particular it allows to distinguish between locally homogeneous (”superstastics” type) and rapidly varying (”diffusing diffusivity” type) environments (see, e.g. figure 2). It also makes it possible to show how the heterogeneity of the original Langevin model induces the non-linear memory structure of the studied process, which is visible in the simulated data.

The work is structured as follows: in section 2 we provide more information on the process of Brownian yet non-Gaussian diffusion along with a primer to autoregressive models. Section 3 then introduces an intuitive physical derivation of the autoregressive model. The situation when the correlation times of the random diffusion coefficients is comparatively short is then considered in detail in section 4. In section 5 we show how the statistical methods of time series can be used to qualify and quantify the non-Gaussianity. There, we also use the analytical formulas for moments to detect the non-linear memory structure specific to the considered model. We summarise our results and discuss their utility in section 6. The paper closes with an appendix, in which the derivation of the moments is presented.

2 Physical stochastic modelling and autoregressive models

2.1 Brownian yet non-Gaussian diffusion

In the original works on Brownian yet non-Gaussian diffusion the linear MSD (1) was observed along with the Laplace shape (4) of the PDF. Other findings suggest the presence of more general stretched Gaussian tails

[TABLE]

with the stretching parameter $\delta\neq 2$. PDFs of this type describe diffusion that can be normal ($\alpha=1$) or anomalous of the form (3) with $\alpha\neq 1$. It was found that $\delta=1$ and $\alpha$ between 0.75 and 0.25 for tracer diffusion in living bacteria and eukaryotic cells [32]. Similarly, in simulations of crowded membranes it was shown that $\delta$ is between 1.3 and 1.6, and $\alpha$ below 1 [33].111Gamma distributions were found in [26], see the modelling in [34].

This behaviour can be explained by the fact that the considered medium is spatially or temporally heterogeneous [23, 24, 28, 35, 36, 37]. In such a system the thermal fluctuations are still Gaussian but the observed displacements are mixtures of these Gaussian contributions with random weights, effecting the non-Gaussian outcome. In practice, this is commonly realised by two classes of diffusion models, namely by so-called "superstatistics" and "diffusing diffusivity" models.

Superstatistics is a term proposed by Cohen and Beck [38, 39, 40] and stands for superposition of statistics. It refers to a model with random parameters, that are fixed for each trajectory. This is a form of hierarchical or multilevel modelling [41], which is also close to the Bayesian inference method. The distributions which arise this way are called compound or mixture [42]. In the context of diffusion this mixture approach can be explained in the following way: each observed trajectory either moves in its own neighbourhood, that has distinct properties affecting the motion. These features are supposed to not vary significantly at temporal and spatial scales specific to the observed trajectories. If the motion is not confined, this assumption must be viewed as a short-time approximation. Another possible origin is the situation when the diffusion characteristics of the observed particles vary from one to another, for instance, the radius of different tracer beads. In this case the compound distributions are observed at long times, as well, but this possibility can be excluded for experiments in which the type of the diffusing object is precisely controlled.

The simplest type of a mixture model is Brownian diffusion with a randomised, but global diffusion coefficient $D$. In it, the infinitesimal increments of the position process are given by

[TABLE]

where $B$ represents standard Brownian motion. If we assume that $D$ has an exponential distribution, the emerging process is characterised by the Laplace PDF (4). Similarly, a stretched exponential PDF of $D$ leads to relation (5).222Of course, other forms may also emerge from such an approach, for instance, Gamma [34] or stable [28] distributions. The randomness of $D$ was confirmed in many experiments [26, 32, 37] and also in numerical studies [33]. Since there is no fundamental reason to pinpoint the diffusion coefficient as the only heterogeneous property of the system, other parameters can also be made superstatistical. For example, a random viscosity or memory kernel in the Langevin equation also lead to non-Gaussian PDFs and generally non-linear MSD [43]. It must be stressed that by introducing a parameter which is random and fixed trajectory-wise, we break not only homogeneity, but also locality and independence. By sharing the same random parameter, the memory structure of the process becomes non-linear and non-Gaussian [43]. This dependence is strong enough to render the resulting motion non-ergodic.

This non-ergodic property is not always desired and can be avoided when the system parameters are allowed to be both random and time-dependent, as in the diffusing diffusivity model. Fixed, deterministic values are replaced by a stochastic variation of the parameter, and for this reason the resulting models are called doubly stochastic [44]. If the introduced parameter evolution is ergodic, the resulting model will be ergodic, but still non-Gaussian. For example, one can consider a time-varying diffusion coefficient leading to the diffusing diffusivity approach proposed by Chubinsky and Slater [45].333In this sense the diffusing diffusivity approach is similar in kind to the diffusing waiting times in correlated continuous time random walks [47, 48]. In this model small Brownian displacements $\mathrm{d}B(t)$ are assumed to be randomly rescaled by random $D(t)$, such that the resulting motion is described by the stochastic equation444Here and below such formulas should be understood as stochastic $L^{2}$ integrals and equations [46].

[TABLE]

Further details of the process depend on the specific choice of $D(t)$—the analysis is often complicated. However, an important case was solved by Cox, Ingerson and Ross [49] who proposed a stochastic differential equation with suitable mean-reverting property that leads to $D(t)$ with exponential type of memory. Their work was motivated by financial applications and was expressed in the language of option prices and volatility modelling. More recently, this approach was applied to physically relevant quantities in [28, 35, 50, 51, 52], in particular, for non-stationary initial conditions [34].

2.2 Autoregressive models

In contrast to models originating from physics, autoregressive models have their roots in filtering, optimal control theory, and economics [29]. In this context the measured sequence of observations $V_{k}$ is assumed to be an output of a system, which is linked linearly to the white Gaussian noise input sequence $Z_{k}$. The classical autoregressive moving average (ARMA) processes [30, 31] are defined by the recursive relation

[TABLE]

This process is denoted by ARMA$(p,q)$ and thus explicitly contains the range parameters $p$ and $q$. The term autoregressive AR($p$) reflects the linear dependence of the observed $V_{k}$ on $p$ past values $V_{k-1},\ldots,V_{k-p}$, and the moving average MA($q$) represents a linear combination of the last value of the noise $Z_{k}$ and $q$ previous values $Z_{k-1},\ldots,Z_{k-q}$. The admissible coefficients $\phi^{i}$ and $\theta^{j}$ are chosen such that the sequence $V_{k}$ is stationary, resembling the physical velocity process (hence the choice of notation "$V$") or highly confined motion.555For non-stationary processes the ARIMA model (”I” stands for ”integrated”) may be used: for ARIMA($p,1,q$) the differences $\Delta V_{k}=V_{k}-V_{k-1}$ are assumed to be ARMA$(p,q)$, for ARIMA($p,n,q$) the $n$th differences are ARMA$(p,q)$.

In most of the literature the model (8) was not meant to explain the behaviour of the system. Rather, the philosophy was concentrated on controlling the data. This can be achieved by finding a reasonably small set of $\phi^{i}$ and $\theta^{j}$ that sufficiently well describe the observed memory structure in the data. The procedure typically employs methods based on mean-square optimisation and information theory [30, 31]. Given the estimates of the $\phi^{i}$ and $\theta^{j}$ the future behaviour can be predicted (in the sense of the least square error or similar), which is a crucial element, for instance, in financial forecasting. Additionally, by inverting relation (8), the $Z_{k}$ can be estimated from the $V_{k}$, therefore various white noise tests can be used to verify the goodness of fit of the model [31].

It is in fact remarkable how many concepts of non-equilibrium and biological physics have their independent counterparts in time series analysis. For example, the study of anomalous diffusion and long-range dependence can be performed by using autoregressive fractionally integrated (ARFIMA) models where the white noise sequence $Z_{k}$ is replaced by the fractional noise $Z_{k}^{d}$ [53, 54] with a power-law covariance structure determined by the parameter $d$, $\mathbb{E}[Z_{k}^{d}Z_{k+j}^{d}]\propto j^{2d-1}$ (one can think of $d$ as Hurst exponent minus 1/2). It can be explicitly defined as the result of the power-law filter applied to $Z_{k}$,

[TABLE]

For ARMA processes, the left hand side of (8)—the AR part— is responsible for modelling an exponential decay of the covariance, whereas the MA part introduces short, finite time corrections. Additional filtering of the white noise presented in (9) extends the ARMA model to the ARFIMA with power-law memory, which was very successfully applied to various anomalous diffusion phenomena encoded in biological data [55]. An analogue of Lévy flights and diffusion with a power-law PDF is the ARMA model with stable noise, where the Gaussian $Z_{k}$ are replaced by non-Gaussian stable random variables [56]. Similarly, fractional Lévy stable motion corresponds to ARFIMA with stable noise [57, 54].

Another line of development in time series analysis, originating from Box-Jenkins models, are non-linear generalisations, mainly the so-called autoregressive conditional heteroscedasticity (ARCH) [58] and generalised ARCH (GARCH) models [59].666A set of random variables is heteroscedastic—as opposed to homoscedastic—if there is at least one sub-population that has different variability from the rest, where ”variability” could be measured in terms of the variance or other dispersion measures. These allow for the parameterisation and prediction of a non-constant variance and proved to be very useful for financial time series; their invention prompted the bestowal of the Nobel Memorial Prize in Economic Sciences in 2003 to Granger and Engle. When integrated ARCH-type processes exhibit non-Gaussian distributions for short times along with a linear time dependence of the MSD. Precisely this observation leads us to pursue the question whether there is a fundamental connection between the two worlds of time series analysis using autoregressive methods and physical models of Brownian yet non-Gaussian type.

The aforementioned notion of heteroscedasticity is related to the time-dependence of the conditional variance [60, 61]. The basic ARCH$(q$) model is defined as

[TABLE]

If an ARMA model is assumed for the noise variance, the model becomes a GARCH$(p,q)$ model,

[TABLE]

A GARCH model has an ARMA-type representation, so that many of its properties are similar to those of ARMA models, for instance, we can estimate the GARCH parameters by the same technique as for ARMA processes [29]. In most applications it is enough to consider the order of the model as $(1,1)$. ARFIMA combined with GARCH describes both power-law decay of the correlation function with finite-lag effects (ARFIMA part) and varying diffusion parameter (GARCH part) [62]. For example, it can describe inhomogeneous diffusion in the cell membrane [63] or solar X-ray variability [64]. While the resemblance to diffusing diffusivity models (7) is indisputable: $V_{k}$ models velocity and $\Sigma_{k}$ models the evolving diffusion coefficient, there is to date no physical derivation of relations (10) and (11).

Another approach related to heteroscedasticity, non-Gaussianity and linear MSD is the main point of interest in this work: we assume that the linear coefficients $\phi^{i}$ and $\theta^{j}$ in the definition of ARMA (8) are replaced by random variables $\Phi_{k}^{i}$ and $\Theta_{k}^{j}$ that are independent of the noise $Z_{k}$:

[TABLE]

This is the class of doubly stochastic models, random coefficient ARMA (rcARMA) [65, 66].

3 Physical derivation of the autoregressive model

Let us start with the classical Langevin equation for the velocity of a diffusing particle,

[TABLE]

which is Newton’s second law with the Stokes dissipative force $-\beta V$ and the stochastic force $F$. In the classical setting $F$ is given by white Gaussian noise with amplitude $k_{B}T\beta$ determined by the fluctuation-dissipation relation [67, 68]. Heterogeneity of the medium can be modelled by making the parameters of equation (13) time-dependent and random—then they describe the local state of the environment of the diffusing particle. All possible correlation effects caused by the particle repeatedly visiting the same areas of the phase space are assumed to be reflected by the memory structure of random functions modelling these local parameters. Starting from the next section we will neglect this dependence entirely, which corresponds to an annealed picture that is also inherent in the current physical diffusing diffusivity models.This assumption can be also justified when the environment changes continuously. The resulting equation can then be written in the form

[TABLE]

We do not provide any more fundamental derivation behind this formula, for us $\Lambda(t)$ and $D(t)$ are effective parameters whose physical meaning is determined by equation (14) alone: $D(t)$ describes the local effective amplitude of velocity gains and $\Lambda(t)$ the local effective linear damping or relaxation of the velocity. In financial language these parameters would be called stochastic return rate and stochastic volatility. Take note this class of models differs from the superstatistical approach of Beck and Cohen [38], who assume the Langevin dynamics and stochasticity of environement can be considered separately. In (30) we consider a particular case leading to superstatistics.

Formula (14) resembles the rcARMA. Indeed, in the standard Euler scheme of numerical simulations, one fixes a time delay $\Delta t$ between observations, denotes $V_{k}=V(k\Delta t)$ and approximates $\mathrm{d}V(t)\approx V_{k}-V_{k-1}$, which results in

[TABLE]

This is an rcAR(1) process with AR coefficient $1-\Lambda_{k}\Delta t$. It is easy to see that the same reasoning applies to any linear stochastic differential equation and thus the classical schemes of discrete-time simulations return rcARMA processes approximating the continuous-time solutions.

In fact, the connections between rcARMA and models of diffusing diffusivity reach deeper than that. Classical results from the theory of time series analysis establish that there is an exact correspondence between ARMA and solutions of linear stochastic equations with constant coefficients—no approximation is needed [69]. The classical Langevin equation for a position of a particle confined in an harmonic potential,

[TABLE]

fits into this category. The resulting discrete motion is the ARMA(2,1) process

[TABLE]

with AR(2) coefficients

[TABLE]

The MA(1) coefficient $\theta$ can be easily calculated numerically or expressed by a somewhat complicated but elementary formula [70]. This class of relations has direct application to the statistical analysis of the data and can be used, for instance, to calculate the exact discrete-time power spectral density (the spectrum given by the Fourier series of $X_{k}$) without the need of the commonly used approximation $\sum_{k}f(X_{k})\Delta t\approx\int f(X(t))\mathrm{d}t$ [71].

The core of the physical interpretation we propose in this work is that the derivation used to obtain (16) and (17) can be generalised for the case of linear equations with time-dependent coefficients. Namely, we multiply (14) by the integrating factor $\exp(\int\Lambda(t)\mathrm{d}t)$ and integrate from $(k-1)\Delta t$ to $k\Delta t$, obtaining

[TABLE]

After some rearrangement it takes the sleek form

[TABLE]

with

[TABLE]

The random coefficient AR form of the left hand side is clearly present, but it is not immediately clear what is the structure of $Z_{k}$ on the right. Different $Z_{k}$ use increments $\mathrm{d}B(s)$ from disjoint intervals $[(k-1)\Delta t,k\Delta t]$. The Gaussian increments $\mathrm{d}B(s)$ are stationary and do not depend on $k$ but they are rescaled by $D(t)$ and the exponent of $\Lambda(t)$, which makes their conditional variances a sequence of random variables,

[TABLE]

Conditioned on $\Lambda$ and $D$ the variables $Z_{k}$ are Gaussian, independent, and have the variance specified above. Thus the $Z_{k}$ can be decomposed as $Z_{k}=\Theta_{k}W_{k}$ for a Gaussian white noise $W_{k}$ and coefficients \Theta_{k}=\sqrt{\mathbb{E}\mathopen{}\mathclose{{}\left[Z_{k}^{2}|\Lambda,D}\right]}, which is a function of random $\Lambda(t)$ and $D(t)$. This is but exactly the rcARMA

[TABLE]

in which the regressive coefficients $0\leq\Phi_{k}\leq 1$ model the local return or relaxation rates, and the amplitude coefficients $\Theta_{k}\geq 0$ describe the heteroscedasticity of fluctuations, that is the variability of the fluctuations’ dispersion. An illustration how the changes of $\Phi_{k}$ are visible in the data can be seen in figure 1.

Each $\Phi_{k}$ is a decreasing functional of the values $\Lambda(t)$ from the interval $[(k-1)\Delta t,k\Delta t]$. Similarly $\Theta_{k}$ is increasing with $D(t)$ and decreasing with $\Lambda(t)$. Therefore $\Phi_{k}$ and $\Theta_{k}$ are two sequences which should mirror the memory structure of $\Lambda(t)$ and $D(t)$. They are dependent, because they both are functionals of the same range of values $\Lambda(t)$. In general the form of this dependence is highly non-linear. When the parameters can be assumed to not vary significantly in the interval $\Delta t$, $\Lambda(t)\approx\Lambda_{k}$ and $D(t)\approx D_{k}$, the approximate relation reads

[TABLE]

In practical applications this formula can be further simplified. Namely, the function $(1-x^{2})/(-\ln(x))$ is approximately linear on the interval $[0,1]$, thus a reasonable approximation is $\Theta_{k}\approx\sqrt{D_{k}\Phi_{k}}$.

No matter what is the distribution of the $\Phi_{k}$ and $\Theta_{k}$, the solution of (23) itself can be expressed in a simple manner. Repeating the recursive relation between $V_{k}$ and $V_{k-1}$ we can write the explicit formula linking $V_{k}$ and $V_{k-n}$ for any $n$ in the form

[TABLE]

Away from the degenerate case $\Lambda(t)\equiv 0$ it is true that $\Phi_{k}<1$, therefore when $n\to\infty$ the series converges towards

[TABLE]

We see that each $V_{k}$ is a mixture of Gaussian variables $W_{k}$ with random weights. The result is not Gaussian, but can be characterised as conditionally Gaussian $\mathcal{N}(0,S^{2})$ with conditional variance

[TABLE]

The distribution of $S^{2}$ determines the non-Gaussianity: the more spread out $S^{2}$ is, the more the PDFs of the velocity and displacement differ from the Gaussian shape. This property can be expressed in terms of the excess kurtosis, which is directly related to the relative standard deviation (RSD) of $S^{2}$. Namely,

[TABLE]

is three times the RSD squared. Different variants of the excess kurtosis or the RSD of $S^{2}$ are sometimes called "non-Gaussianity parameter" [72, 73], but here we will use more precise terminology. Because both are positive, the distribution of $V$ is always leptokurtic, that is, it has tails thicker than a Gaussian. It is worth stressing that the kurtosis is only one of many possible measures of non-Gaussianity, but it is a convenient one because of the easy estimation.

The distribution of $S^{2}$ is, in general, hard to analyse even for simple models of $\Phi$ and $\Theta$. One exception is encountered when the coefficients stay constant for long intervals of time $T_{\alpha}$ (precisely, \mathbb{E}\mathopen{}\mathclose{{}\left[T_{\alpha}}\right] larger than the mean relaxation time 1/\mathbb{E}\mathopen{}\mathclose{{}\left[\Phi}\right]). In such a situation $S^{2}$ spends short times transitioning between these periods and therefore in the majority they attain an equilibrium value obtained from taking the coefficients in (27) constant,

[TABLE]

where $\alpha$ is an identification of the interval. The resulting distribution of $S^{2}$ is thus

[TABLE]

with $\widetilde{\Theta}$ and $\widetilde{\Phi}$ being random variables constructed from the values $\Theta_{\alpha},\Phi_{\alpha},\Theta_{\beta},\Phi_{\beta},\Theta_{\gamma},$ $\Phi_{\gamma},\ldots$ with probabilities weighted by the relative mean lengths of the intervals $T_{\alpha},T_{\beta},T_{\gamma},...$. In this approximation one trajectory is effectively the same as an ensemble of equal-length trajectories, each with its own $\widetilde{\Theta},\widetilde{\Phi}$, which are glued together at the ends, merging into a single trajectory. (This is essentially a form of superstatistics.) For a simple demonstration see figure 2. There one of the limiting cases are long periods of constant coefficients, as in (30), but one can also observe the second limit of rapidly varying coefficients. This is the subject of the next section.

4 Short memory random coefficient AR(1)

When the correlation time of $\Lambda(t)$ and $D(t)$ is much smaller than the observation time, $t_{c}\ll\Delta t$, the integrals defining $\Phi_{k},\Phi_{k+1}$ and $\Theta_{k},\Theta_{k+1}$ contain only the ratio $t_{c}/\Delta t$ of the total mass corresponding to dependent values of $\Lambda(t),D(t)$. Thus, $\Phi_{k}$ and $\Theta_{k}$ can be assumed to be sequences of independent and identically distributed (i.i.d.) variables. Note that for fixed $k$ the pair $\Phi_{k},\Theta_{k}$ may, and generally will be dependent. Averaging (27) yields

[TABLE]

Similarly, \mathbb{E}\mathopen{}\mathclose{{}\left[S^{4}}\right] can be calculated, as detailed in the appendix. The moments \mathbb{E}\mathopen{}\mathclose{{}\left[S^{2}}\right] and \mathbb{E}\mathopen{}\mathclose{{}\left[S^{4}}\right] together determine the kurtosis of $V$. Accordingly, the motion is again non-Gaussian. The covariance also has a simple form: from the recursive formula (3) we can easily derive the geometric decay

[TABLE]

The form of the covariance also determines the MSD: for the velocity it converges exponentially to a constant $\delta^{2}_{V}(j)=2(1-r_{V}(j))$ and for the position $X_{k}\mathrel{\mathop{:}}=\int_{0}^{k\Delta t}V(\tau)\mathrm{d}\tau$ it is linear for $j\gg 1$,

[TABLE]

A similar formula can be obtained if one tries to recreate displacements given $V_{k}$: then $\widetilde{X}_{k}\mathrel{\mathop{:}}=\Delta t\sum_{i=1}^{k}V_{i}$ leads to the formula

[TABLE]

A small error in determining the slope of the true $X(t)$ is made in this procedure, but the general behaviour of the MSD is the same. Note that the shapes of all the above functions depend only on the mean value of $\Phi_{k}$, this is a diffusion with covariance and MSD (up to a scaling factor) as if the system were homogenised—however, at the same time it is non-Gaussian, see figure 3.

As a concrete example of a physical model which leads to this type of process let us assume that the diffusivity is constant (with no loss of generality we take $D(t)=1$). Then we consider a diffusing particle interacting with high viscosity traps distributed in a Poissonian manner, that is, the waiting time between subsequent trapping events have the exponential distribution $\mathcal{E}(\rho)$. By a high-viscosity trap we mean a short-time interval when $\Lambda(t)=\infty$, which immediately kills the momentum of the particle—this is a form of stochastic resetting [74]. Outside of these events the damping rate is constant $\Lambda(t)=\lambda$. The probability that the particle does not meet any traps between times $(k-1)\Delta t$ and $k\Delta t$ is $\exp(-\rho\Delta t)$: in this case $\Phi_{k}=\exp(-\lambda\Delta t)$ and $\Theta_{k}=\sqrt{(1-\exp(-2\lambda\Delta t))/(2\lambda)}$. The probability that the particle meets at least one trap in an interval $[(k-1)\Delta t,k\Delta t]$ is of probability $1-\exp(-\rho\Delta t)$, and the corresponding AR(1) coefficient $\Phi_{k}$ becomes 0. The discrete noise $Z_{k}$ for this $k$ has a more complex structure. All Gaussian fluctuations before the last trapping event are erased, which is visible from (21). However the remaining fluctuations still count. Denoting by $T_{k}$ the difference between the last trapping event and $k\Delta t$ we obtain the formula for $\Theta_{k}^{\prime}$ corresponding to the event $\Phi_{k}=0$,

[TABLE]

This corresponds to a small correction $\mathcal{O}(\Delta t)$, which makes the resulting PDF smooth (otherwise the trapping events would be visible as a Dirac delta at $x=0$). The approximation on the right holds for $\lambda\Delta t\ll 1$. The relation between $\Lambda(t)$ and $\Phi_{k},\Theta_{k}$ is shown in figure 4 using an exemplary trajectory.

Variables $T_{k}$ which determine $\Theta_{k}$ have the distribution of an exponential variable conditioned by the requirement that it has value lower than $\Delta t$ (that is, the trapping event occurred in the considered interval). It has the PDF

[TABLE]

As we determined the exact distribution of $\Phi_{k}$ and $\Theta_{k}$, the formula for the covariance (32) can be made completely explicit, namely

[TABLE]

The variance \mathbb{E}\mathopen{}\mathclose{{}\left[S^{2}}\right] has a complicated form, but the geometric decay rate is given by $(\rho+\lambda)\Delta t$. It is an intuitive result: meeting a trap erases the momentum of the particle, which thus loses all connection to its history at a frequency proportional to the density of traps.

Studying the distribution of $V_{k}$ in general is hard, even for i.i.d. coefficients, because the conditional variance (27) is given by an infinite sum with terms which are dependent, as the same $\Phi_{k-i}$ reappear in them. In the studied example the situation is specific and simpler, because the sum is actually finite. The coefficients $\Phi_{k-i}$ are a series of Bernoulli trials—after a series of $j-1$ non-zero ones appearing with probability $\exp(-\rho\Delta t)^{j}$ the first zero occurs with probability $1-\exp(-\rho\Delta t)$ and the summation stops. Conditioned by this event the first $j-1$ values $\Theta_{k-i}$ are deterministic, but the last term in the sum, $\Theta_{k-j}^{\prime}$ has the conditional distribution (35). This leads to the formula for characteristic function,

[TABLE]

Here the Laplace transform of the conditional $\Theta_{k-j}^{\prime}\overset{d}{=}\Theta_{k-j}|\{\Phi_{k-j}=0\}$ is given by

[TABLE]

The last approximation works for $\lambda\Delta t\ll 1$ which follows from (35). The limiting case $\lambda\Delta t\approx 0$ is also interesting: then $\Phi_{k}$ are i.i.d. variables with Bernoulli distribution $\mathrm{P}(\Phi_{k}=1)=\exp(-\rho\Delta t)=1-\mathrm{P}(\Phi_{k}=0)$, as usual $\Theta_{k}=\sqrt{\Phi_{k}}$—for simplicity we may neglect the fluctuations counted in the next interval after the trapping events. Directly from $(\ref{eq:S2})$ we see that $S^{2}$ has a geometric distribution. Any such variable can be expressed as an integer part of some exponentially distributed variable, $E$ say, and therefore it has values between $E$ and $E-1$. As mentioned in the introduction, the mixture of Gaussian variables with variance $E$ has exponential tails [75], so this is the case for $S^{2}$, as well. This also provides an intuitive argument for the presence of the tails with thickness between Gaussian and Laplace in the general case considered in this section (see figure 3). As the $\Phi_{k}$ are bounded by $1$, the tails of $S^{2}$ correspond to long stretches of $\Phi_{k}$ close to $1$, and as they are i.i.d. the probability of such an event decays geometrically with the number of $\Phi_{k}$.

5 Statistical analysis

We now address the statistical procedures in the context of stochastic processes with random coefficients.

5.1 Memory

In (32) we showed that the covariance of i.i.d. rcAR(1) is, up to a scale factor, the same as for the AR(1) with coefficient \phi=\mathbb{E}\mathopen{}\mathclose{{}\left[\Phi_{k}}\right]. The same is also true for general i.i.d. rcARMA, it can be demonstrated by multiplying both sides of the equation by past values of the process and deriving the recursive formulas for the covariance (these are called Yule-Walker equations)—they differ from the non-random coefficient case only by the single equation which determines the variance [76]. On the one hand it makes the additional randomness harder to detect, but on the other hand it allows to use powerful ARMA based methods of analysis. As we mentioned in section 2.2, AR (thus also rcAR) have a covariance function that is a mixture of exponentials, and the MA part is responsible for short time corrections. For the continuous-time processes the only option of statistical verification would seem to be fitting the estimated covariance function and checking the validity of the result. For rcARMA a more intuitive tool is available: the partial autocorrelation function. Given a stationary time series $X_{k}$ it measures the linear dependence between $X_{k}$ and $X_{k-j}$ when an in-between set of variables $\mathcal{I}_{k,j}\mathrel{\mathop{:}}=X_{k-1},X_{k_{2}},\ldots X_{k-j+1}$ is removed. Denoting the least-squares projection operator by $\mathcal{P}$ it can be defined as

[TABLE]

The partial autocorrelation can be effectively estimated using standard methods implemented in the popular statistical packages (such as Python Statsmodels, R tseries, Julia StatsBase, or Matlab Econometrics Toolbox)—the usual command is the standard abbreviation pacf. It is intuitively clear that $\alpha_{X}$ measures the "direct" dependence between $X_{k}$ and $X_{k-j}$ which for AR(p) and rcAR(p), given their recursive definition, should be non-zero for the first $p$ values and zero further on. For the MA part one can reverse the definition and express the noise $Z_{k}$ as the geometrical sum of $X_{k}$, hence, the partial autocorrelation decays geometrically [30, 31]. For the full ARMA model the two effects are additive.

For the i.i.d. rcAR(1) the statistical procedure becomes very simple: just estimate $\alpha_{X}$ and check if it has a significant non-zero value only for $j=1$. This value is simply \alpha_{X}(1)=\mathbb{E}\mathopen{}\mathclose{{}\left[\Phi_{k}}\right]. For an illustration see figure 5. For higher order AR processes there are also exact equations linking the values of partial autocorrelation, covariance and coefficients of the model [30, 31].

We now proceed to show that there actually is a way to detect differences between the dependence structure of ARMA and rcARMA. First, we estimate the mean rcAR(1) coefficient \varphi=\mathbb{E}\mathopen{}\mathclose{{}\left[\Phi_{k}}\right], for instance, by using the partial autocorrelation. Next, the typical method of statistical verification for AR processes is to consider

[TABLE]

For the AR(1) with deterministic coefficient $\phi=\varphi$ this would be an estimate of the noise $Z_{k}$. Therefore, one could then use many of the testing methods if the series is i.i.d. However, for the rcAR(1),

[TABLE]

As \mathbb{E}\mathopen{}\mathclose{{}\left[\Phi_{k}-\varphi}\right]=0 the series is still uncorrelated, $\mathbb{E}[\widetilde{Z}_{k}\widetilde{Z}_{k+1}]=0$. It is a natural consequence of the applied procedure, which fits the system by removing the linear dependence. But the series $\widetilde{Z}_{k}$ remains non-linearly dependent. The value $\widetilde{Z}_{k}$ still depends on $V_{k-1}$ and consequently also on $\Phi_{k-1}$ appearing in the expression for $\widetilde{Z}_{k-1}$. This purely non-linear dependence can be detected only using non-linear measures of memory. The codifference function (see [77] for a detailed discussion) is one possible choice. It was proposed as a tool to study $\alpha$-stable variables, because it is finite even for variables with no second, or even lower, moments [78]—however its usefulness is by no means limited to this class. One of the few possible close definitions is the formula

[TABLE]

For any Gaussian process the codifference equals the covariance, so the difference between the two indicate any non-Gaussianity of the PDF and, what is important here, of the memory structure. For i.i.d. rcAR it can be expressed as a function of coefficients, and it can be analytically approximated, see (51) and the derivation below. As a simple example let us take the uniformly distributed $\Phi_{k}\overset{d}{=}\mathcal{U}(0,1)$ and $\Theta_{k}=\sqrt{D_{k}\Phi_{k}},D_{k}\overset{d}{=}\mathcal{U}(0,1)$. The covariance is zero as $\mathbb{E}[\widetilde{Z}_{k}\widetilde{Z}_{+1}]=0$, but the simulation yields $\tau_{\widetilde{Z}}^{1}(1)\approx 0.0023$, in good agreement with the analytic approximation is close $412/170535=0.0024$. Two more examples of rcAR(1) with higher dispersion of i.i.d. coefficients and even more prominent non-linear memory are shown in figure 6.

5.2 Non-Gaussianity

One possible approach of visualising the non-Gaussian behaviour of a diffusive process is to quantify the non-Gaussian nature of its dispersion in a manner similar to how we quantified non-Gaussian memory in the last section. We propose to use the non-linear measure of dispersion [77]

[TABLE]

which is called log-characteristic function (LCF). It has the property that if the process is Gaussian it equals the MSD for any $\theta$. Therefore estimating it and plotting it together with the MSD reveals the non-Gaussian nature of the observed process. As the MSD averages over squares of the data, small values become even smaller and larger values become even more emphasised. Thus the MSD gives a large importance to the spread of extremal values. Conversely, for the LCF large values are translated into rapid oscillations of the cosine function and cancel each other, the LCF measures mainly the spread of the bulk.

For the classical case of non-Gaussian normal diffusion (4) the LCF is linear for small times, but logarithmic for large times, so it gets dominated by the linear MSD. This insight is more general. Random coefficient models discussed in this work can only describe diffusion processes with tails thicker than Gaussian. This means that compared to Gaussian diffusion with the same MSD more mass is located in the tails than in the bulk, and the LCF must be smaller than the MSD.

To use the LCF in practice one needs to fix the parameter $\theta$. For too small values the LCF converges to the MSD because $\cos(\theta x)\approx 1-\frac{1}{2}\theta^{2}x^{2}$, for too large values the estimation becomes challenging, because $\ln(0)$ diverges. In a typical case reasonable values are those for which \theta\mathbb{E}\mathopen{}\mathclose{{}\left[|X(t)|}\right]\approx\pi so that $X(t)$ explores mainly the first period of the cosine function. In figure 7 we show a practical application of this technique performed using a relatively small sample of 500 integrated rcAR trajectories.

For rigorous statistical testing of the Gaussianity we propose to use the standard Jarque-Bera (JB) test [79]. It is a goodness of fit test of whether sample data have skewness and kurtosis that match a normal distribution. For the sample $\{x_{1},\ldots,x_{n}\}$ of observation the JB statistic is defined as

[TABLE]

where $\bar{s}$ is the sample skewness and $\bar{\kappa}$ is the sample excess kurtosis. Samples from a normal distribution have an expected skewness of [math] and an expected kurtosis of $3$ (excess kurtosis 0). Any deviation from this increases the JB statistic. The test is considered as standard and is implemented in various numerical packages, such as R tseries, Python SciPy.stats, Julia HypothesisTests, or Matlab.

This test is considered as one of the most powerful tests on normality but, if non-Gaussianity is detected, it does provide information on the origins of this behaviour. However, the idea of observing the kurtosis leads to an algorithm to distinguish between Gaussian and non-Gaussian (in particular, Lévy stable) distributions [80]. It is based on the empirical cumulative excess kurtosis (ECEK) which is defined as follows,

[TABLE]

where $\bar{x}$ is an arithmetic mean of the random sample. This simple statistic serves as an indicator for whether there is a noticeable difference between Gaussian and non-Gaussian distributions. For the Gaussian case, for large numbers of observations it converges to the theoretical excess kurtosis 0, while for the non-Gaussian case the ECEK does not tend to 0 with increasing number $k$ of observations and, moreover, for distributions that do not have a finite forth moment it does not converge at all and behaves chaotically [80].

To distinguish between Gaussian and non-Gaussian distributions we also advocate the application of the discrimination algorithm based on examining the rate of convergence to the Gaussian law by means of the central limit theorem [81]. The idea of the algorithm is to analyse the convergence of the estimated index $\alpha$ of stability for sequential bootstrapped samples from the analysed data. If the estimated values converge to 2, then the data are light-tailed and belong to the domain of attraction of the Gaussian law. In particular, if the data are Gaussian, the estimated values should be equal to 2 for most of the cases. If the data belong to the domain of attraction of a non-Gaussian Lévy stable law, then the values should converge to a constant less than 2.

To illustrate the usefulness of the above statistical tools we consider an rcAR process for which the $\Phi_{k}$ are independent uniformly distributed on the interval $(0,0.95)$ and $\Theta_{k}=\sqrt{\Phi_{k}}$. In the left panel of figure 8 we can see a simulated trajectory of $V_{k}$ of length 1000 points. In the middle panel we display a plot of the ECEK for the simulated trajectory. We can see that the function converges to a constant close to 2 which clearly indicates a non-Gaussian behaviour and a finite positive excess kurtosis of the underlying distribution (a leptokurtic distribution). We also calculated analytically the exact value of the excess kurtosis (equations (28) and (50)) and obtained the value $1.83$, which coincides with the final values of the ECEK function. Finally, in the right panel of figure 8 we show the estimated values of the stability index $\alpha$, for different non-overlapping consecutive blocks of length $M$. For the first sample ($M=1$; corresponding to the whole trajectory), the value is significantly lower than 2, but for the other $M$ the values increase and tend to 2. This suggests that the simulated trajectory exhibits a non-Gaussian behaviour but its distribution belongs to the normal domain of attraction. The obtained results clearly suggest that the simulated trajectory exhibits a non-Gaussian behaviour and the underlying distribution is heavier-tailed than Gaussian (is leptokurtic) but belongs to the normal domain of attraction. The non-Gaussianity is also confirmed by the JB test which returns a $p$-value less than $0.001$.

5.3 Parameter estimation for rcAR processes

Finding reliable estimators of the parameters of the rcARMA models is a very important challenge. In this section we present preliminary ideas for some special cases of the models.

For an rcAR process with $\Phi$ being a simple function one can estimate its parameters either by extracting the constant periods of $\Phi$ and applying classical fitting techniques for the autoregressive processes to the extracted parts, or to consider a method which directly incorporates the information on switching between different autoregressive processes.

We now follow the latter idea and apply an algorithm based on hidden Markov Models (HMMs) [64, 82, 83]. This algorithm assumes that the trajectories switch among discrete diffusive states according to a stochastic (Markov) process. We consider a modification of the classical HMM, where the Markov regime switching is combined with AR$(1)$ processes [84]. To show the usefulness of this technique we take into account a two-regime parameter switching model, that is, a model with both regimes driven by AR$(1)$ processes with the autoregressive parameters $\phi_{1}$ and $\phi_{2}$. In figure 9 we present a simulated trajectory of $V_{k}$ of length 1000 for which $\Phi_{\alpha}$ changes from 0.2 to 0.8 at point 400. The estimated regime switching point is $399$ which is very close to the true value, and the estimated $\phi_{1}$ and $\phi_{2}$ coefficients are $0.17$, $0.79$ which confirms the usefulness of the procedure.

6 Discussion

As modern experiments such as single particle tracking routinely produce large amounts of data for the thermally and actively driven motion of test particles, the extraction of physical information from the garnered data becomes ever more important. On the one hand, this is met by the analysis of a growing number of complementary observables such as time averaged MSD [85, 86], higher order moments or mean maximal statistics [87], p-variation methods [88], first-passage methods [89], or single trajectory power spectra [8, 90]. On the other hand objective methods such as maximum likelihood approaches [91, 92, 93, 94] or machine learning [95] are being recognised as useful tools. Here the goal is to identify the physical model giving rise to the observed motion and estimate the involved parameters.

An alternative method, well established in financial market studies, is time series analysis. The latter seemingly does not share a lot with the physical approach of modelling. It puts the emphasis not on explaining the observed systems, but instead concentrates on fitting and prediction. Given a class of discrete linear filters suited to describe the general features of the data (for us the non-Gaussianity and linear MSD) it offers a wide choice of powerful methods for finding the optimal model and validating it. In other words by design it concentrates on the phenomenological description.

While a pure time series analysis approach initially may be discouraging for physically-minded researchers, the undeniable effectiveness of the time series methods make any connections to the framework of physics interesting and important. As noted in section 3 the correspondence between the Gaussian Langevin equation and the ARMA process is known (in fact, it was established as early as in 1959 [96]), but it never seemed to have gained a widespread appeal in the physical or biological physics community. Nevertheless, some recent progress was made, especially in topics related to the identification of fractional dynamics [63, 97] and accounting for the distortions made by measurement devices [71, 98].

This paper continues this line of research and aims at promoting to sample the best of both worlds. This is hard to achieve on a general level, but we study here important cases when a fruitful compromise can be made, and new information provided in the challenging analysis of non-Gaussian diffusion. Omitting the dynamics at time scales significantly shorter than what is available from observation (shorter than the sampling time $\Delta t$) leads to conceptually simple autoregressive models. These are straightforward and quick to simulate (using a medium class computer we were generating ten million values in under 2 seconds) and they provide a wide choice of statistical methods. Some of those, such as the partial autocorrelation, are intrinsically linked to discrete dynamics and have no continuous-time counterpart. Others, such as the kurtosis and codifference, profit from the analytic methods available for discrete variables.

In section 5 we illustrated important properties of the rcARMA processes and provided a list of statistical tools that can be useful in their identification and validation. In particular we showed that the codifference function can be used to distinguish between ARMA and rcARMA processes. To illustrate the non-Gaussianity of the rcARMA models we considered a non-linear measure of dispersion, namely the log-characteristic function and the empirical cumulative excess kurtosis, and studied the domain of attraction of the underlying distribution. Finally, we presented a hidden Markov model methodology which can be useful in fitting the process parameters. We also stress that the presented tools can be successfully applied for quite moderate lengths of the trajectories.

Our considerations are general and not limited to any particular system, as we discuss for related models based on random coefficient autoregression, which stem from a widely used form of the Langevin equation and reasonable models of heterogeneous environment. Thus, the provided methods are also universal. We hope that this study will promote the use of autoregressive models in modern data analysis, as well as prompt further studies into the physical meaning of these models.

JŚ acknowledges support through the Polish National Science Centre grant 2016/22/M/ST1/00233. KB acknowledges support through Beethoven Grant No. DFG-NCN 2016/23/G/ST1/04083. RM acknowledges funding from DFG, grants ME1535/6-1 and ME1535/7-1. RM was supported by an Alexander von Humboldt Polish Honorary Research Scholarship from the Foundation for Polish Science (Fundacja na rzecz Nauki Polski).

Appendix A Appendix: Moments of $S$ and the approximation of codifference

Here we show the calculation of \mathbb{E}\mathopen{}\mathclose{{}\left[S^{4}}\right] and the approximation of the codifference for an rcAR with i.i.d. coefficients. Starting with (27) we write

[TABLE]

We express this square as the sum of the elements squared plus twice the sum of each element multiplied by all elements to the right of it. The sum of squares is simply \mathbb{E}\mathopen{}\mathclose{{}\left[\Theta_{k}^{4}}\right]/(1-\mathbb{E}\mathopen{}\mathclose{{}\left[\Phi_{k}^{4}}\right]). As for the rest, in each element a factor of type \mathbb{E}\mathopen{}\mathclose{{}\left[\Phi_{k-i}^{2}\Theta_{k-i}^{2}}\right]\mathbb{E}\mathopen{}\mathclose{{}\left[\Theta_{k-j}^{2}}\right]=\mathbb{E}\mathopen{}\mathclose{{}\left[\Phi_{k}^{2}\Theta_{k}^{2}}\right]\mathbb{E}\mathopen{}\mathclose{{}\left[\Theta_{k}^{2}}\right] appears. Note that it does not necessarily decouple, because in general $\Phi_{k}$ and $\Theta_{k}$ for fixed $k$ may be dependent. The remaining factor is

[TABLE]

We recognise the product of two geometric series in the above. Finally, the fourth moment becomes

[TABLE]

Given these moments, one can easily calculate the kurtosis. It also makes it possible to obtain an analytic approximation of the codifference of the estimated noise \widetilde{Z}_{k}\mathrel{\mathop{:}}=V_{k}-\varphi V_{k-1},\varphi=\mathbb{E}\mathopen{}\mathclose{{}\left[\Phi_{k}}\right]. We show the procedure for $\tau_{\widetilde{Z}}^{\theta}(1)$. The method for other $\tau_{\widetilde{Z}}^{\theta}(j),j>1$ would be the same. We denote $\widetilde{\Phi}_{k}\mathrel{\mathop{:}}=\Phi_{k}-\varphi$ and start with a straightforward conditioning which yields

[TABLE]

where $S^{2}$ is assumed to have the distribution as in (27) and be independent from $\Phi_{k},\Theta_{k},\Phi_{k-1},\Theta_{k-1}$. Expanding the exponents in a Taylor series yields an approximation for the codifference which depends only on the moments of $\Phi_{k}$ and $\Theta_{k}$. Expansion up to first order, $\rme^{-x}\approx 1-x$, is equivalent to calculating the covariance, so it is zero. The second order, $\rme^{-x}\approx 1-x+x^{2}/2$ in general already provides quite good non-zero estimate. For the numerator it is

[TABLE]

The corresponding denominator is expanded into

[TABLE]

The expected values that appear in this formulas are just linear combinations of the first four moments of $\Phi_{k},\Theta_{k}$ and their products. For variables generated as powers of a uniform distribution these are just integrals over polynomials. Also, without significant change the logarithm can be replaced by the fraction $\ln(x)\approx x-1$, which leads to a formula being a fraction of coefficient moments up to the fourth order.