Nowadays, the concept of a Dynamical System covers systems of any nature (physical, chemical, biological, economical etc) both deterministic and stochastic. Such a system can be described with differential equations, functions from algebra of logic, graphs, Markov chains, etc. (Butenin et al., 1987, p.8)
Let us consider an ensamble of N different phase trajectories (starting from the same initial condition). This ensamble is determined by an ensamble of realizations of random sources perturbing the system. The statistical ensamble of N → ∞ realizations defines a stochastic process {X(t), t ∈ T}. A single phase trajectory (a realization of stochastic process) is denoted as x(t).
Let us measure realizations of the process X(t): x_{1}, x_{2}, … x_{n} at times t_{1}, t_{2}, … t_{n} to find the joint probability density p_{n}(x_{1}, t_{1}; …; x_{n}, t_{n}) for their occurence. The complete information on the dynamics of the stochastic process X(t) is given by the infinite sequence of ntime (ndimensional) joint probability densities:
where δ() is the Dirac's deltafunction. Since an infinite number of distributions law is needed to describe fully a random process it is impossible to do this, in general. Thus one should extract as much knowledge as possible from an one  dimensional distribution and loworder moments (expectation, variance and autocorrelation function).
A numerical approximation of a trajectory x(t) is given by such a finite set of random variables
that these variables x^{τ}_{i} approximate unknown random variables x(t_{i}) at
where
A numerical approximation of a stochastic process X(t) is given by such a set of K random variables (set of trajectories)
Generally stationarity means invariance of some property under a time shift.
A stochastic process is called strictly stationary (stationarity in a narrow sense) if the joint probability density depends on the differences between two subsequent times. A change of t_{i} → t_{i} + τ for all i does not affect the probability density. Neither characteristic of a process changes under a time shift.
A stochastic process is called weakly stationary (stationarity in a wide sense) if its expectation, variance and autocorrelation function do not change under a time shift.
If a property of a stochastic process (e.g. its variance) does not change in time then the process is stationary with regard to that property.
An asymptotically stationary stochastic process is when a change of t_{i} → t_{i} + τ for t_{i} → ∞ does not affect the probability density.
Note. For a stationary process (in both senses), expectation and variance are constant values
and autocorrelation function depends only on the time difference
pquantile line (ζ_{p}) is a deterministic line used to estimate dynamic behavior of time series. For p from the range (0, 1) and a given stochastic process X
a pquantile line is defined as
and is approximated by the following sequence obtained from K realization of the stochastic process
where x_{i,(ki(p))} is a random value taken from a non  decaying sequence built from K realization of a stochastic process at the ith time
k_{i}(p) =  { 

Note. For a stationary processes p  quantile line are parallel to each other.
A process is called strictly ergodic if all its characteristics can be determined from its single (infinitely long) realization. Temporal averaging and ensamble averaging give the same result. In the case when only certain characteristics of a process can be restored from the single realization, the process is called ergodic in respect to those characteristics. An ergodic process is stationary. The reverse statement is not true, in general.
Spectral analysis means any representation of a signal as a superposition of some
Fourier Transform means decomposition of a signal into harmonic components. A detected signal may represent a superposition of signals from several sources. When each of those sources demonstrates harmonic oscillations with its own frequency then its intensity in observed signal is reflected by the value of the power spectrum at the corresponding frequency.
The Fourier Transform of a continuous signal is defined as
∞  
F(f) =  ∫  x(t) e^{i2πft} dt 
∞ 
In the case of discrete time a random sequence {x(t_{n})}^{N1}_{n=0} is converting into frequency sequence via the discrete Fourier transform
N1  
DFT(m) =  ∑  x(n)e^{i2πnm/N}. 
n=0 
DFT(m) is the mth component of DFT, N is a number of input signal time points and x(n) is the nth time point of time series x. The DFT sequence is of complex type and its absolute value reads
DFT is computed using the standard algorithm of the Fast Fourier Transform (FFT) with a smoothing window:
N1  
DFT_{w}(m) =  ∑  w(n)x(n)e^{i2πnm/N} 
n=0 
where w(n) window can be of
triangular type
w(n) =  { 

Hanning type
Hamming type
The Power Spectrum (PS) of a single deterministic periodic signal is welldefined. However, when one considers a realisation of a random process X(t) the problem is more complicated. First of all, random quantities are almost always nonperiodic. Moreover, integrals for continuous Fourier transforms almost always do not exist. Then spectral properties of a stochastic process are described by the finitary Fourier Transform i. e. integrals over an interval [T/2, T/2]:
T/2  
FT_{T}(ω) =  ∫  x(t) e^{(iωt)} dt 
T/2 
To find an estimator for the power spectrum S(ω) the expectation of the FT_{T}(ω)^{2} is computed thus the power spectrum is given by the expression
S(ω) =  lim  2 ⟨ FT_{T}(ω) ^{2}⟩ / T 
T → ∞ 
In the program, the set of {DFT_{abs}(m)}^{M}_{m=1} values is obtained numerically with the discrete Fourier transform (DFT) for each of N realizations. Then to get an estimator of the power spectrum S(m) averaging over realizations of the considered process is performed.
Note 1. The Power Spectrum allow to detect different sources of oscillations and estimate their relative intensity.
Note 2. Fourier transforms and Power Spectrum of stochastic processes exist only as expressions averaged values over different realizations.
Note 3. The power spectrum of stationary processes equals the doubled Fourier transform of the autocorrelation function (the WienerKhinchin theorem).
The firstorder moment is the expectation of X
∞  
E[X(t)] =  ∫  x p(x,t) dx 
∞ 
where p(x,t) is a probability density function.
An unbiased estimator of the expectation from a sample of N independent values is the sample mean:
N  
⟨ x ⟩ = 1/N  ∑  x_{i} 
i=1 
Thus, the expectation E[X(t_{m})] for fixed t_{m} is computed as arithmetic means over N realizations.
The secondorder central moment is called variance
∞  
E[(X(t)  E[X(t)])^{2}] =  ∫  (x  E[X(t)])^{2} p(x,t) dx 
∞ 
where p(x,t) is a probability density function.
An unbiased estimator of the variance from a sample of N independent values is the sample variance:
N  
σ^{2} = 1/(N  1)  ∑  (x_{i}  ⟨x⟩)^{2} 
i=1 
Hence, the expectation E[(X(t_{m})  E[X(t_{m})])_{2}] for fixed t_{m} is computed as arithmetic means over N realizations.
The correlation function of a stochastic process X(t) can be computed in two ways. First, in both cases, N realizations of the signal X(t) are detected or loaded from a file. Then
for fixed t_{1} and t_{2} the expectations are computed as arithmetic means over N realizations
N  
E[X(t_{1,2})] = 1/N  ∑  x_{k}(t_{1,2}) 
k=1 
After that, ensemble averaging is done to find the autocovariance function defined as the expectation
If it is normalized by rootmeansquared deviations one gets autocorrelation function
where
N  
σ^{2}(t_{1,2}) = E[(X(t_{1,2})  E[X(t_{1,2})])^{2}] = 1/(N1)  ∑  (x_{k}(t_{1,2})  E[X(t_{1,2})])^{2} 
k=1 
the timeaveraged value of the x variable for a single realization of an ergodic process X(t) is computed as an arithmetic mean
I1  
⟨x⟩_{t} = 1/I  ∑  x(t_{i}) 
i=0 
where i is a moment in time. After that, time averaging is done to find the mean product of
Im1  
ρ(τ) = ⟨x(t)x(t+τ)⟩_{t} = 1/(Im)  ∑  x_{i} x_{i+m} 
i=0 
where m = 0,1,…, (I  L), k = 1, …, N and Δt is a time step. L (cutoff) is set as I/2 but can be changed by an user, however, must be greater than 1. And normalized to the maximum value at τ = 0
The cross  correlation function of two stochastic processes X(t) and Y(t) is computed to reveal character of coupling between sources of both signals.
To do this, first, N realizations of signals X(t) and Y(t) are detected or loaded from a file. Next,
ensemble averaging is done to find the crosscovariance function defined as the expectation
If it is normalized by rootmeansquared deviations one gets crosscorrelation function
the timeaveraged value of x and y variables for a single realization of ergodic processes X(t) and Y(t) are computed
I1  
⟨s⟩_{t} = 1/I  ∑  s(t_{i}) 
i=0 
where i denotes temporal index and s = {x, y}. After that, time averaging is done to find the mean product of
Im1  
ρ^{xy}(τ) = ⟨x(t)y(t+τ)⟩_{t} ≈ 1/(Im)  ∑  x_{i} y_{i+m}  = ρ^{xy}_{m} 
i=0 
and normalized to the maximum value
Note 1. The absolute value of a correlation function between two processes vary between 0 and 1. For deterministically linear dependence between them the correlation function reaches 1 while for statistically independent processes it is 0.
Note 2. For the stationary processes the absolute value of the correlation function obeys the CauchySchwartz inequality
where
ρ^{xy}(τ) =  ∫  dxdyxyp_{xy}(x, t; y, t + τ). 
The correlation time τ_{C} for a continuous normalized autocorrelation function ρ_{N}(τ) is defined as
∞  
τ_{C} =  ∫  dτρ_{N}(τ) 
0 
whereas for a discrete sequence {ρ_{N,m}} one has
M  
τ_{C} =  ∑  Δtρ_{N,m} 
m=0 
where Δt is a time step.
Fourier transform of
the autocorrelation function:
∞  
G_{xx}(ω) =  ∫  dτ ρ^{xx}(τ)e^{iωτ} 
∞ 
and of the crosscorrelation function
∞  
G_{xy}(ω) =  ∫  dτ ρ^{xy}(τ)e^{iωτ} 
∞ 
exists if the correlation functions are absolutely integrable.
Probability Density Function (PDF) of a random probe of K size embeded in the range [a, b] can be estimated by H(x) function defined as
K  
H_{K}(x) = 1/(Kh)  ∑  I_{(xl, xl+1]}(X_{k}) 
k=1 
for x ∈ (x_{l}, x_{l+1}] where l = 0,1,…, L1 and
Probability Density Function (PDF) can be estimated by f_{K}(x) function defined as
K  
f_{K}(t_{i},x_{j}) = 1/K  ∑  J((x_{j}  X^{τ}_{i,k})/b_{Ki})/b_{Ki } 
k=1 
where x_{j} ∈ ℜ and t_{i} ∈ [t_{0}, T]. J denotes a RosenblattParzen's kernel of
rectangular type
triangular type
gaussian type
optimal type
and
is a random sample of K size from an unknown distribution. τ is a time step and time indices i form the sequence i = 0,…,I.
Entropy function H(t) is defined as
∞  
H(t) =   ∫  f(t,x) log f(t,x) dx 
∞ 
where f(t,x) is a probability density function. H(t) is approximated by the sum
J  
H_{i} = H(t_{i}) =   ∑  f_{i,j} log f_{i,j} Δx 
j=0 
Stochastic exponent X(t) of a semimartingale process Z(t) is defined as the solution of the equation
t  
X(t) = 1 +  ∫  X(s) dZ(s) 
0 
where Z(s) can be
gaussian process with the random measure
where N(0,1) is the standard normal distribution
α  stable Levy motion with the random measure
where S_{α, β} is a random variable from αstable Levy distribution
poissonian process with the random measure
where N_{λ, i} is a compensated poissonian process with intensity λ
or their linear combination
Z(s) often is a cadlag process (or a RCLL process i. e. Right Continuous with Left Limits) i.e.
X(t) process where t ∈ [0, T] is approximated by a discrete process
defined as
where X^{τ}_{0} = 1, i = 0, …,I, and M([t_{i1},t_{i})) is a random measure defined above for each stochastic process.
Let's define set of points
and corresponding ranges
then the expression
defines probability that values of the X process belong to the given range.
Package for machine learning  OptFinderML.