Numerical Simulation of the MRW Process

The bottleneck of numerical simulation of the MRW process (8) is simulation of logarithmically correlated noise! Дг [k]. Simulation of the discrete Gaussian noise process with given autocorrelation function (covariance matrix) is a well – known problem and is subjected to a trade-off: exact simulation processes usually requires a lot of computation resources, and fast algorithms typically provide only approximated solution. The most known exact simulation method is based on the Cholesky or LU-decomposition of the covariance matrix into lower – and upper-triangle matrices (Davis 1987). Though this method is very efficient for short time-series, having computational complexity of O(N2) it is not suitable for simulation of long (e. g. N > 103) realizations. Much more efficient decomposition algorithm is the so-called Circulant Embedding Method (CEM) (Dietrich and Newsam 1997) that is based on the Fast Fourier Transform (FFT) and thus has complexity of O(N log N).

Consider N x N covariance matrix R with elements Rp>q = r [|p — q|] = r [k] for k = 1,…,N — 1, where r[k] is the required covariance matrix and in our case is given by (10). CEM consist in embedding of matrix R into a larger 2M x 2M matrix S (where M > N — 1). The optimal case of M = N — 1 is called the minimal embedding. The first row of S, denoted by s, consists of two parts of length N — 1 each: the original first row of R following with the first row of R in reverse order:

Подпись: s [k] = r [k], s [2M — k] = r [k], Подпись: k = 0 ,...,N — 1, k = 1,...,N — 2.
Подпись: (15)

Since matrix S is circulant, any matrix extracted along its diagonal is a copy of R. One of the properties of circulant matrices tells that matrix S can be decomposed into a product S = FDFH, where F is a matrix of discrete Fourier transform coefficients, and FH is the conjugate transpose of F. The matrix D is diagonal, and elements along the diagonal can be obtained by the discrete Fourier transform of first row or column of S: s = Fs.

Let us construct a vector y = FD1/2x, where x is a complex vector, having iid N(0,1/ real and imaginary parts. It was shown (Dietrich and Newsam 1997) that any vector of size N, extracted from either real or imaginary part of y has covariance matrix exactly equal to R. It should be noted that necessary and sufficient condition of existence of auxiliary vector y is the matrix S to be nonnegative definite (nonnegative eigenvalues sm of the matrix S). Though for strictly positive definite Toeplix matrices R the existence of nonnegative circulant embedding was proven (Dembo et al. 1989), there is no general recipe and the necessary and sufficient condition should be tested for any specific covariance matrix R.

Подпись: sm Подпись: N-1 rk exp k=0 Подпись: mk 2л i . 2 (N — 1) Подпись: 2 (N-1)-1 + ^2 r2(N-1 )-k exp k=N image166

Without providing a general proof that circulant matrix S for covariance matrix given by (10) is nonnegative definite, while performing simulations we have tested numerically eigenvalues S for all used parameter sets. According to the property of circulant matrices, eigenvalues sm of matrix S are equal to the discrete Fourier transform of the first row of S:

n=X[10] log Д+(— 1)m log

Подпись: s image168 Подпись: (16)

where i[11] 2 [12] [13] = —1 and rk = Cov [юд-[і],!Д-[і + k]] are given by (10). After the rearrangement and substitution of (10) we obtain:

As one can see, X2 is a multiplicative coefficient in (16) and the sign of sm is fully determined by the relation between LjД – and N. For any used combination of L/Д – and N we have tested that values (16) are nonnegative: sm > 0.

The CEM algorithm for minimal embedding case could be summarized in the following steps (Dietrich and Newsam 1997):

5. Generate complex iid N(0,1) vector x of length 2(N — 1);

6. Multiply result from step 4 by x;

7. Evaluate y via inverse FFT of the result of previous step;

8. Extract vector of length N from real part of y;

9. Another vector could be extracted from real part of y;

10. To generate additional realizations, proceed to step 5.

The computational complexity of this algorithm is O(N log N) compared to O(N2) for the Cholesky-based method. Additionally, the memory requirements is O(N) for CEM compared to O(N2) for Cholesky decomposition. This allows to generate extremely long (N = 1020) realizations of MRW random process using standard Core i5 4Gb RAM computer.

Leave a reply

You may use these HTML tags and attributes: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <s> <strike> <strong>