Skip to content

val-link/LongMemoryProcess

Folders and files

NameName
Last commit message
Last commit date

Latest commit

 

History

8 Commits
 
 
 
 
 
 

Repository files navigation

Long Memory Process

This repo contais a Julia script to generate realizations of a class of stationary Gaussian stochastic processes with extremely long memory. In detail the script can be used to generate realizations X ( t ) of a real stationary Gaussian process with correlation function

cf ( t ) = E ( X ( t 0 + t ) X ( t 0 ) ) = ( Γ ( y + 1 , 0 ) Γ ( y + 1 , θ max | t | ) ) ( θ max | t | ) y 1

with θ max > 0 and y > 1 . Γ is the incomplete gamma function. The correlation function decays algebraically for large t (heavy tailed)

cf ( t ) Γ ( y + 1 ) ( θ max | t | ) y 1 .

Theory

Standard tools to generate stationary processes can become slow when the memory time is very large. The scheme that I use here approximates the true process as a finite sum of Ornstein-Uhlenbeck processes. This may be much more efficient than convolution in the time domain (as in ARFIMA(0,d,0)) or using a Fourier method. The advantage of the specific form for the correlation function is that it has the following integral representation

cf ( t ) = ( θ max ) y 1 0 θ max θ y e θ | t | d θ .

θ max determines the fastest decay rate that we allow for and can be used to modify the short time behavior.

When we apply a numerical quadrature scheme to the integral it reduces to a finite sum over exponentials (in t )

cf ( t ) ( θ max ) y 1 k w k ( θ k ) y e θ k | t |

with weights w k and abscissas θ k . Exponentially decaying correlation functions are realized by a stationary Ornstein-Uhlenbeck process for which the transition probability is known analytically. Thus, the numerical quadrature scheme delivers an approximation of the process as a sum over a finite number of stationary Ornstein-Uhlenbeck processes. An efficient scheme for the diverging integrand is tanh-sinh quadrature (this is used in the code).

This strategy allows for an efficient generation of realizations of stationary Gaussian processes with very long memory (i.e. y close to 1 ).

Usage

Generation of realizations

First create a LongMemoryProcess struct.

lmp = LongMemoryProcess(y, θ_max)

To generate a single realization along a vector of times t one can then simply call

Xₜ = lmp(t)

Accuracy

You can check the accuracy by comparing the true correlation function with the fitted function (that is the correlation function of the approximated process).

y = -0.9; θ_max = 1.0; t_f = 1e8;
lmp = LongMemoryProcess(y, θ_max);
relative_error = abs(cf(y, θ_max, t_f) - cf_fit(lmp, t_f)) / cf(y, θ_max, t_f)   # 9.013557696260914e-7

Calling cf requires the incomplete gamma function via SpecialFunctions.jl. If you need more speed and less accuracy you can decrease the number of Ornstein-Uhlenbeck processes with the optional argument n (the default value is n = 401).

number_of_OU_processes = 61;
lmp = LongMemoryProcess(y, θ_max, n = number_of_OU_processes);
relative_error = abs(cf(y, θ_max, t_f) - cf_fit(lmp, t_f)) / cf(y, θ_max, t_f)   # 0.0035874732157112664

Advanced stuff

If you want to plug in the process as noise in a differential equation you can use an SDE solver and integrate the Ornstein-Uhlenbeck processes along. The individual Ornstein-Uhlenbeck processes obey the SDE

d x k ( t ) = θ k x k ( t ) d t + σ k d W k ( t ) ,

and the full process is just their sum X ( t ) = k x k ( t ) . The parameters σ k , θ k are fields of the LongMemoryProcess struct.

θ = lmp.θ; σ = lmp.σ;

To obtain a valid initial condition for x k ( t ) simply generate a realization for a single time step.

X₀, x₀ = lmp(Vector([0.]); return_full=true);
X₀[end] == sum(x₀)   # true

Citing

If you found this project useful for your work please directly cite this GitHub repository.

About

Efficiently generate Gaussian stochastic processes with heavy-tailed algebraic correlations.

Topics

Resources

License

Stars

Watchers

Forks

Languages