Simulation of circularly-symmetric fractional Gaussian noises and modulated fractional Gaussian noises Jean-François Coeurjolly and Emilio Porcu May 2017 The file simGCP.R contains four different functions related to the manuscript: • • • •

find.mt: lists highly compposite numbers larger than 2n − 1 for a given n. simCFGN: simulation of a circularly-symmetric fractional Gaussian noise. simModFGN: simulation of the sum of modulated fractional Gaussian noises. efficiencyPlot: image plot in terms of parameter values of the validity of the assumption of nonnegativity of the circulant matrix C.

It also contains two other R functions which are devoted to estimate the Hurst exponent of a circularlysymmetric fractional Gaussian noise. These two functions (dilation, estCFBM) are not commented hereafter.

The function find.mt We recommend to embed the covariance matrix into a circulant matrix C with size (m, e m) e where m e is a highly composite number of the form 3i 5j 7k 11l where i, j, k, l are integers larger than 2n − 1. Such numbers up to max.mt=1e7 (default value) are obtained as follows > vec.mt=find.mt(1000) > head(vec.mt) ## ## ## ## ## ## ##

[1,] [2,] [3,] [4,] [5,] [6,]

i 4 3 7 2 0 2

j 2 0 0 1 0 2

k 0 1 0 2 4 0

l mt=3^i 5^j 7^k 11^l 0 2025 1 2079 0 2187 0 2205 0 2401 1 2475

Hence, the first highly composite number larger than 2n − 1 = 1999 is 2025.

Simulation of a circularly symmetric Gaussian process To generate a sample path with length n = 1000, Hurst exponent H = .7 and η parameter η = 2/3 tan(πH) (remind that the condition for which the model is well-defined is |η| ≤ | tan(πH)|)), with circulant matrix C with dimension (m, e m) e with m e = 1999 (m e must be larger than 2n − 1), type > + > > > >

out=simCFGN(n=1000,H=.8,eta=2/3*tan(pi*.8),mt=1999, allowTruncation=FALSE) ## plot of the cumulated sums par(mfrow=c(1,2)) plot(cumsum(Re(out$Z)),type="l",ylab="Real part");grid() plot(cumsum(Im(out$Z)),type="l",ylab="Imaginary part");grid()

1

250

0

0

50

100

Imaginary part

150

200

200 150 100 50

Real part

0

200

600

1000

0

Index

200

600

1000

Index

> par(mfrow=c(1,1)) If you want the first highly composite number larger than 2n − 1, type mt=NULL which is the default. > out=simCFGN(n=1000,H=.8,eta=.3,mt=NULL, + allowTruncation=FALSE) ## No entry for mt...use the first highly composite integer larger than 2n-1, mt= 2025 If you try with H = 0.95 and η = 3/4 tan(πH) for instance and keep allowTruncation=FALSE this will produce an error since some eigenvalues are negative and ask you to allow the truncation of eigenvalues. The approximate simulation is described in Appendix B in the manuscript. A message gives the number of eigenvalues which are ≤ 0. A plot of the bound of P (k∆k∞ > x) is provided if the parameter plot.error is set to TRUE. > out=simCFGN(n=1000,H=.95,eta=3/4*tan(pi*0.95),mt=NULL, + allowTruncation=TRUE) ## No entry for mt...use the first highly composite integer larger than 2n-1, mt= 2025 ## 373 / 2025 eigenvalues are negative...use of the approximate simulation > out=simCFGN(n=1000,H=.95,eta=3/4*tan(pi*0.95),mt=NULL, + allowTruncation=TRUE,plot.error=TRUE) ## No entry for mt...use the first highly composite integer larger than 2n-1, mt= 2025 ## 373 / 2025 eigenvalues are negative...use of the approximate simulation

2

1.0

Bound of the supremum norm of the error

0.6 0.4 0.0

0.2

P( ∆ ∞ > x) ≤ e(x)

0.8

e−1(0.1)=0.9 e−1(0.05)=0.94 e−1(0.01)=1.01

0.7

0.8

0.9

1.0

1.1

x The output of the function simCFGN is a list with the following objects: • • • • • •

Z: the realization (or approximate realization) of the CFGN; complex vector. lc.before: vector of length m e of the eigenvalues of C before truncation. lC: vector of truncated eigenvalues. delta: parameter δ2 defined in Appendix B.1. nb.neg: number of eigenvalues which are negative. r: vector containing the values γ(0), . . . , γ(n − 1) where γ is the covariance function of the circularlysymmetric fractional Gaussian noise.

A plot of the eigenvalues is thus easily produced.

Simulation of the sum of modulated fGn The function simModFGN is quite similar to the function simFGN in terms of paramters and outputs. The parameter H and eta are replaced by vH and vphi where vH and vphi are respectively the vector of Hurst exponents and modulations. > out= simModFGN(n=1000,vH=.4,vphi=pi/2,mt=NULL, + allowTruncation=TRUE,plot.error=TRUE) ## No entry for mt...use the first highly composite integer larger than 2n-1, mt= 2025 Such parameters lead to an exact simulation which is not the case for the following choice: > out= simModFGN(n=1000,vH=.95,vphi=pi/2,mt=NULL, + allowTruncation=TRUE,plot.error=TRUE) ## No entry for mt...use the first highly composite integer larger than 2n-1, mt= 2025 ## 1012 / 2025 eigenvalues are negative...use of the approximate simulation

3

1.0

Bound of the supremum norm of the error

0.6 0.4 0.0

0.2

P( ∆ ∞ > x) ≤ e(x)

0.8

e−1(0.1)=2.23 e−1(0.05)=2.32 e−1(0.01)=2.52

1.6

1.8

2.0

2.2

2.4

2.6

2.8

x The sum of two modulated fGn is done as follows: > out= simModFGN(n=1000,vH=c(.4,.7),vphi=c(pi/4,pi/2),mt=NULL, + allowTruncation=TRUE,plot.error=TRUE) ## No entry for mt...use the first highly composite integer larger than 2n-1, mt= 2025

Efficiency plots The function efficiencyPlot allows to investigate parameter values (H, η) for a circularly-symmetric fGn and (H, φ) for a modulated fGn for which the circulant matrix C is non-negative. The main parameter is mt: • If mt
4

> out=efficiencyPlot(n=1000,mt=NULL,choice='CFGN',l=200,rangeH=c(.51,.99), + rangePar=c(-1,1)*abs(tan(pi*.51)),verbose=FALSE,posLeg=c(.85,-.3))

30

~ =2025 n=1000, m

10

0.20

η

0.15

−30

−10

δ=1 δ=1 δ=1 δ=2 δ=3 δ=1

0.6

0.7

0.8

4 3 2 3 4

0.10 0.05 0.00

0.9

H (blue surface= 80.7 %; >0 eigenv.) The blue surface means that all eigenvalues were larger or equal to zero for the first highly composite number larger than 2n − 1. As underlined by the legend, 80.7% of the parameter values yield that situation. Let us now zoom on the interval [.75, .99]: > out=efficiencyPlot(n=1000,mt=NULL,choice='CFGN',l=200,rangeH=c(.75,.99), + rangePar=c(-1,1)*abs(tan(pi*.75)),verbose=FALSE,posLeg=c(.92,0))

1.0

~ =2025 n=1000, m

0.0

0.20 0.15 δ=1 δ=1 δ=1 δ=2 δ=3 δ=1

−1.0 −0.5

η

0.5

0.25

0.75

0.80

0.85

0.90

4 3 2 3 4

0.10 0.05 0.00

0.95

H (blue surface= 69.5 %; >0 eigenv.) We can observe that when η ≤ 2/3| tan(πH)| there are no negative eigenvalues. When negative values are

5

observed, we notice that for a large set of (H, η) the fraction of negative values is lower than 5%. Let us now have a look at what happens when we increase m: e > out=efficiencyPlot(n=1000,mt=3^9,choice='CFGN',l=200,rangeH=c(.75,.99), + rangePar=c(-1,1)*abs(tan(pi*.75)),verbose=FALSE,posLeg=c(.92,0))

1.0

~ =19683 n=1000, m

0.0

0.20 0.15 δ=1 δ=1 δ=1 δ=2 δ=3 δ=1

−1.0 −0.5

η

0.5

0.25

0.75

0.80

0.85

0.90

4 3 2 3 4

0.10 0.05 0.00

0.95

H (blue surface= 69.5 %; >0 eigenv.) We can see that there is no real effect on the “blue set” of parameter values. The gain is that the other admissible values seem to provide a fraction of negative values which tends to decrease with m. e The following instruction investigates a (single) modulated fractional Gaussian noise. > out=efficiencyPlot(n=1000,mt=NULL,choice='ModFGN',l=300,rangeH=c(.01,.99), + rangePar=c(0,pi),verbose=FALSE)

6

3.0

~ =2025 n=1000, m

2.0

0.4 φ

0.3

1.0

0.2 0.1

0.0

0.0 0.2

0.4

0.6

0.8

H (blue surface= 89.5 %; >0 eigenv.) Let us zoom on the interval [.8, .99] for H and on a smaller interval for φ. > out=efficiencyPlot(n=1000,mt=NULL,choice='ModFGN',l=200,rangeH=c(.8,.99), + rangePar=c(-.004,.004),verbose=FALSE)

0.004

~ =2025 n=1000, m

0.3

0.000

φ

0.4

0.2

−0.004

0.1 0.0 0.80

0.85

0.90

0.95

H (blue surface= 49.2 %; >0 eigenv.) Here is what happens when m e is larger. > out=efficiencyPlot(n=1000,mt=3^8,choice='ModFGN',l=200,rangeH=c(.8,.99), + rangePar=c(-.004,.004),verbose=FALSE)

7

0.004

~ =6561 n=1000, m

0.3

0.000

φ

0.4

0.2

−0.004

0.1 0.0 0.80

0.85

0.90

0.95

H (blue surface= 47.4 %; >0 eigenv.) Let us now investigate the first highly composite number larger than 2n − 1 and smaller than mt.max such that C is non-negative, i.e. such that all eigenvalues are non-negative. This is done by supplying the parameter mt=Inf. A white pixel inside the admissible set of parameter values means that no m e ≤ mt.max were found. > out=efficiencyPlot(n=1000,mt=Inf,choice='CFGN',l=100,rangeH=c(.75,.99), + rangePar=c(-1,1),verbose=FALSE,mt.max=20000,posLeg=c(.92,0))

0.0

15000 δ=1 δ=1 δ=1 δ=2 δ=3 δ=1

−1.0 −0.5

η

0.5

1.0

~ =Inf (<20000) n=1000, m

0.75

0.80

0.85

0.90

4 3 2 3 4

10000 5000

0.95

H (blue surface= 69.5 %; corresponds to the minimal embedding) From the previous plot, we clearly see what was guessed previously. The first higlhy composite number larger than 2n − 1 ensures the non-negativity of a large set of parameter values and that increasing m e does not help. The next image plot considers the modulated fractional Gaussian noise. Unlike the CFGN model, we observe

8

that an increase of m e can ensure the non-negativity of C. > out=efficiencyPlot(n=1000,mt=Inf,choice='ModFGN',l=100,rangeH=c(.8,.99), + rangePar=c(-.004,.004),verbose=FALSE,mt.max=20000)

0.004

~ =Inf (<20000) n=1000, m

0.000

φ

15000 10000

−0.004

5000

0.80

0.85

0.90

0.95

H (blue surface= 50.3 %; corresponds to the minimal embedding)

Efficiency plots generated for the manuscript

1.0

> out=efficiencyPlot(n=1000,mt=NULL,choice='CFGN',l=200,rangeH=c(.75,.99), + rangePar=c(-1,1)*abs(tan(pi*.75)),verbose=FALSE,posLeg=c(.935,-0.17),paper=TRUE)

0.5

0.25

0.0

0.15

−0.5

δ=1 δ=1 δ=1 δ=2 δ=3 δ=1

−1.0

η

0.20

0.75

0.80

0.85

0.90 H

9

0.95

4 3 2 3 4

0.10 0.05 0.00

0.0

15000

10000 δ=1 δ=1 δ=1 δ=2 δ=3 δ=1

−1.0

−0.5

η

0.5

1.0

> out=efficiencyPlot(n=1000,mt=Inf,choice='CFGN',l=100,rangeH=c(.75,.99), + rangePar=c(-1,1),verbose=FALSE,mt.max=20000,posLeg=c(.935,-.17),paper=TRUE)

0.75

0.80

0.85

0.90

4 3 2 3 4

5000

0.95

H

0.4 0.3

0.000

φ

0.002

0.004

> out=efficiencyPlot(n=1000,mt=NULL,choice='ModFGN',l=200,rangeH=c(.8,.99), + rangePar=c(-.004,.004),verbose=FALSE,paper=TRUE)

0.2

−0.004

0.1 0.0 0.80

0.85

0.90

0.95

H > out=efficiencyPlot(n=1000,mt=Inf,choice='ModFGN',l=100,rangeH=c(.8,.99), + rangePar=c(-.004,.004),verbose=FALSE,mt.max=20000,paper=TRUE)

10

0.004 0.002 0.000

10000

5000 −0.004

φ

15000

0.80

0.85

0.90

0.95

H

11