VMSTA Modern Stochastics: Theory and Applications 2351-6054 2351-6046 2351-6046 VTeXMokslininkų g. 2A, 08412 Vilnius, Lithuania VMSTA59 10.15559/16-VMSTA59 Research Article Simulation paradoxes related to a fractional Brownian motion with small Hurst index MakoginVitaliimakoginv@ukr.net Department of Probability Theory, Statistics and Actuarial Mathematics, Taras Shevchenko National University of Kyiv, 64, Volodymyrska St., 01601 Kyiv, Ukraine 2016 47201632181190 3152016 1962016 2062016 © 2016 The Author(s). Published by VTeX2016 Open access article under the CC BY license.

We consider the simulation of sample paths of a fractional Brownian motion with small values of the Hurst index and estimate the behavior of the expected maximum. We prove that, for each fixed N, the error of approximation Emaxt[0,1]BH(t)Emaxi=1,NBH(i/N) grows rapidly to as the Hurst index tends to 0.

Fractional Brownian motion Monte Carlo simulations expected maximum discrete approximation 65C50 60G22
Introduction

A fractional Brownian motion {BH(t),t0} is a centered Gaussian stochastic process with covariance function E[BH(t)BH(u)]=12(t2H+u2H|tu|2H),t,u0, where H(0,1) is the Hurst index. The fractional Brownian motion is a self-similar process with index H, that is, for any a>0 0$]]>, {BH(t),t0}=d{aHBH(at),t0}, where =d means the equality of finite-dimensional distributions. Due to self-similarity, we have that, for all T>0 0$]]>, {BH(t),t[0,T]}=d{THBH(T1t),t[0,T]}={THBH(s),s[0,1]}. Based on such a invariance of distributions, it is appropriate to investigate the properties of the fractional Brownian motion only over the time interval [0,1] .

In this paper, we consider the behavior of the maximum functional maxt[0,1]BH(t) with small values of Hurst index.

It should be noted that the fractional Brownian motion process with H=1/2 is the Wiener process {W(t),t0} . The distribution of maxt[0,1]W(t) is known. Namely, P(maxt[0,1]W(t)x)=2π0xey2/2dy,x0, and, therefore, E[maxt[0,1]W(t)]=2π.

Many papers are devoted to the distribution of the maximum functional of the fractional Brownian motion, where usually asymptotic properties for large values of time horizon T are considered. For example, Molchan  has found an asymptotic behavior of small-ball probabilities for the maximum of the fractional Brownian motion. Talagrand  obtained lower bounds for the expected maximum of the fractional Brownian motion. In several works, the distribution of the maximum is investigated when the Hurst index H is close to 1/2 . In particular, this case was considered by Sinai  and recently by Delorme and Weise .

Currently, an analytical expression for the distribution of the maximum of the functional Brownian motion remains unknown. Moreover, the exact value of the expectation of such a functional is unknown too.

From the paper of Borovkov et al.  we know the following bounds: 12Hπeln2Emaxt[0,1]BH(t)<16.3H. On the other hand, we may get an approximate value of the expected maximum using Monte Carlo simulations. That is, for sufficiently large N, Emaxt[0,1]BH(t)Emaxi=1,NBH(i/N). The authors of  obtain an upper bound for the error ΔN of approximation (2). Namely, for N21/H , 0ΔN:=Emaxt[0,1]BH(t)Emaxi=1,NBH(i/N) 2lnNNH(1+4NH+0.0074(lnN)3/2).

The implementation of approximation (2) has technical limitations. Due to modern computer capabilities, we assume that N220106 . Under such conditions, inequality (4) is true when H0.05 , and ΔN<11.18 .

In this article, we make Monte Carlo simulations and estimate Emaxi=1,NBH(i/N) . Also, we investigate the behavior of ΔN with small values of the Hurst index H and show that, for a fixed N, the approximation error ΔN+ as H0 . For the rate of this convergence, when N=220 , we prove the inequality ΔN>c1H1/2c2 c_{1}{H}^{-1/2}-c_{2}$]]>, H(0,1) , where the constants c1=0.2055 and c2=3.4452 are calculated numerically. Thus, when the values of H are small, approximation (2) is not appropriate for evaluation of Emaxt[0,1]BH(t) . The article is organized as follows. The first section presents the methodology of computing. The second section presents the results of computing of the expected maximum of the fractional Brownian motion. In the third section, we obtain a lower bound for the error ΔN and calculate the constants c1 and c2 . Methods of approximate calculations Simulation of a vector <inline-formula id="j_vmsta59_ineq_031"><alternatives> <mml:math><mml:msub><mml:mrow><mml:mo mathvariant="normal" fence="true" stretchy="false">(</mml:mo><mml:msup><mml:mrow><mml:mi mathvariant="italic">B</mml:mi></mml:mrow><mml:mrow><mml:mi mathvariant="italic">H</mml:mi></mml:mrow></mml:msup><mml:mo mathvariant="normal" fence="true" stretchy="false">(</mml:mo><mml:mi mathvariant="italic">i</mml:mi><mml:mo mathvariant="normal" stretchy="false">/</mml:mo><mml:mi mathvariant="italic">N</mml:mi><mml:mo mathvariant="normal" fence="true" stretchy="false">)</mml:mo><mml:mo mathvariant="normal" fence="true" stretchy="false">)</mml:mo></mml:mrow><mml:mrow><mml:mn>1</mml:mn><mml:mo stretchy="false">≤</mml:mo><mml:mi mathvariant="italic">i</mml:mi><mml:mo stretchy="false">≤</mml:mo><mml:mi mathvariant="italic">N</mml:mi></mml:mrow></mml:msub></mml:math> <tex-math><![CDATA[$({B}^{H}(i/N))_{1\le i\le N}$]]></tex-math></alternatives></inline-formula> Let us consider briefly the method proposed by Wood and Chan . Let G be the autocovariance matrix of (BH(1/N),,BH(N/N)) . Embed G in a circulant m×m matrix C given by C=c0c1cm1cm1c0cm2c1c2c0, where cj=1N2H(|j1|2H2j2H+(j+1)2H),0jm2,1N2H((mj1)2H2(mj)2H+(mj+1)2H),m2<jm1. Let m=21+ν , where 2ν is the minimum power of 2 not less than N. Then the matrix C allows a representation C=QJQT , where J is a diagonal matrix of eigenvalues of the matrix C, and Q is the unitary matrix with elements (Q)j,k=1mcjexp(2iπjkm),j,k=0,m1. The eigenvalues λk of the matrix C are equal to λk= j=0m1exp(2iπjkm),k=0,m1. Since Q is unitary, we can set Y=QJ1/2QTZ , where ZN(0,Im) . Therefore, we get YN(0,C) . Thus, the distributions of the vectors (Y0,Y0+Y1,,Y0++YN1) and (BH(1/N),,BH(N/N)) coincide. The method of Wood and Chan is exact and has complexity O(Nlog(N)) . A more detailed description of the algorithm, a comparison with other methods of simulation of the fractional Brownian motion, and a program code are contained in the paper . For reasons of optimization of calculations, simulations in the present paper are made by the method of Wood and Chan. The estimate of the mean value Emaxi=1,NBH(i/N) is a sample mean over the sample of size n. That is why the total complexity of the algorithm is O(nNlog(N)) . Clark’s method Instead of generating samples and computing sample means, there exists a method of Clark  for approximating the expected maximum. Due to this method, the first four moments of the random variable max{ξ1,,ξN} , where (ξ1,,ξN) is a Gaussian vector, are calculated approximately. Since the fractional Brownian motion is a Gaussian process, we put (ξ1,,ξN)=(BH(1/N),,BH(N/N)) and apply Clark’s method for approximate computing of Emaxi=1,NBH(i/N) . Let us illustrate the basic idea of Clark’s method of calculating Emax{ξ,η,τ} , where ξ,η,τ are Gaussian distributed. Let ξ,η,τ be Gaussian random variables. Put a=Var(ξ)+Var(η)Cov(ξ,η) and let a>0 0$]]>. Denote α:=(EξEη)/a . Then we have Emax{ξ,η}=Φ(α)Eξ+Φ(α)Eη+aφ(α);E(max{ξ,η})2=Φ(α)Eξ2+Φ(α)Eη2+aφ(α)(Eξ+Eη), where φ(x)=12πexp(x22) and Φ(x)=xϕ(t)dt .

So, the exact value of Emax{ξ,η} is obtained from the previous proposition.

Let ξ,η,τ be Gaussian random variables. Let Corr(τ,ξ) and Corr(τ,η) be known. Then Corr(τ,max{ξ,η})=Var(ξ)Corr(τ,ξ)Φ(α)+Var(η)Corr(τ,η)Φ(α)E(max{ξ,η})2(Emax{ξ,η})2.

For approximate computing Emax{ξ,η,τ}=Emax{τ,max{ξ,η}} , we assume that max{ξ,η} has a Gaussian distribution. In fact, this is not true, but it allows us to apply formula (5) for random variables τ and max{ξ,η} . Thus, iteratively, we can calculate the approximate mean for any finite number of Gaussian random variables.

Computing the expected maximum

In this section, we present results of approximate computing Emaxi=1,NBH(i/N) by generating random samples and applying Clark’s method. Also, we compare the computational results obtained by these two methods.

The values of the Hurst index are taken from the set {104(1+4i),i=0,24}{102i,i=1,9} . The values of N are chosen from the set {2j,j=8,19} . The values of (BH(1/N),BH(2/N),,BH(N/N)) are simulated by the method of Wood and Chan for each pair N,H with the sample size n=1000 . For each element in the sample, we calculate the following functionals: maxi=1,NBH(i/N), 1N i=1NBH(i/N).

Approximation error of <inline-formula id="j_vmsta59_ineq_071"><alternatives> <mml:math><mml:mstyle displaystyle="false"><mml:mfrac><mml:mrow><mml:mn>1</mml:mn></mml:mrow><mml:mrow><mml:mi mathvariant="italic">N</mml:mi></mml:mrow></mml:mfrac></mml:mstyle><mml:msubsup><mml:mrow><mml:mo largeop="false" movablelimits="false">∑</mml:mo></mml:mrow><mml:mrow><mml:mi mathvariant="italic">i</mml:mi><mml:mo>=</mml:mo><mml:mn>0</mml:mn></mml:mrow><mml:mrow><mml:mi mathvariant="italic">N</mml:mi></mml:mrow></mml:msubsup><mml:msup><mml:mrow><mml:mi mathvariant="italic">B</mml:mi></mml:mrow><mml:mrow><mml:mi mathvariant="italic">H</mml:mi></mml:mrow></mml:msup><mml:mo mathvariant="normal" fence="true" stretchy="false">(</mml:mo><mml:mi mathvariant="italic">i</mml:mi><mml:mo mathvariant="normal" stretchy="false">/</mml:mo><mml:mi mathvariant="italic">N</mml:mi><mml:mo mathvariant="normal" fence="true" stretchy="false">)</mml:mo></mml:math> <tex-math><![CDATA[$\frac{1}{N}{\sum _{i=0}^{N}}{B}^{H}(i/N)$]]></tex-math></alternatives></inline-formula>

We compute the sample mean and variance of (7). The values of theoretical moments of (7) are known: E1N i=0NBH(i/N)=0, E(1N i=0NBH(i/N))2=1N2H+2 i=1Ni2H+11(2H+2),N.

Sample means of 1Ni=0NBH(i/N)

Sample variances of 1Ni=0NBH(i/N)

The sample moments of (7) when H={104(1+4i),i=0,24} and N={2j,j=8,19} are presented in Figs. 1 and 2. In the figures, the lines indicate the theoretical moments and confidence intervals corresponding to the reliability of 95%. The data confirm the correctness of calculations of (7) with the reliability of 95% even for small values of H.

Computing functional <inline-formula id="j_vmsta59_ineq_076"><alternatives> <mml:math><mml:msub><mml:mrow><mml:mo movablelimits="false">max</mml:mo></mml:mrow><mml:mrow><mml:mi mathvariant="italic">i</mml:mi><mml:mo>=</mml:mo><mml:mover accent="false"><mml:mrow><mml:mn>1</mml:mn><mml:mo mathvariant="normal">,</mml:mo><mml:mi mathvariant="italic">N</mml:mi></mml:mrow><mml:mo accent="true">‾</mml:mo></mml:mover></mml:mrow></mml:msub><mml:msup><mml:mrow><mml:mi mathvariant="italic">B</mml:mi></mml:mrow><mml:mrow><mml:mi mathvariant="italic">H</mml:mi></mml:mrow></mml:msup><mml:mo mathvariant="normal" fence="true" stretchy="false">(</mml:mo><mml:mi mathvariant="italic">i</mml:mi><mml:mo mathvariant="normal" stretchy="false">/</mml:mo><mml:mi mathvariant="italic">N</mml:mi><mml:mo mathvariant="normal" fence="true" stretchy="false">)</mml:mo></mml:math> <tex-math><![CDATA[$\max _{i=\overline{1,N}}{B}^{H}(i/N)$]]></tex-math></alternatives></inline-formula>

For each pair N,H , we obtain the sample of values of the maximum functional (6) with sample size 1000. For some values of H, the sample means and approximate values of the expected maximum, obtained by Clark’s method, are presented in Table 1.

The approximate values of the expected maximum

 Sample means of (6) Values due to Clark’s method N\H 0.0900 0.0100 0.0013 0.0001 0.0900 0.0100 0.0013 0.0001 28 1.7017 2.0019 1.9897 1.9769 1.1738 1.8691 1.9696 1.9839 29 1.7693 2.0875 2.1602 2.1360 1.1903 1.9991 2.1194 2.1366 210 1.9487 2.2504 2.3047 2.2854 1.1971 2.1193 2.2604 2.2806 211 2.0138 2.4203 2.4446 2.4184 1.1966 2.2310 2.3939 2.4174 212 2.0886 2.5086 2.5948 2.5334 1.1910 2.3351 2.5208 2.5476 213 2.1938 2.6396 2.6934 2.6885 1.1822 2.4327 2.6420 2.6723 214 2.2591 2.7612 2.7829 2.7940 1.1714 2.5242 2.7579 2.7919 215 2.3327 2.8837 2.9452 2.9258 1.1586 2.6104 2.8693 2.9070 216 2.4050 2.9973 3.0526 3.0464 1.1436 2.6917 2.9765 3.0181 217 2.4620 3.0791 3.1386 3.1121 1.1263 2.7685 3.0798 3.1256 218 2.5328 3.1900 3.2102 3.2421 1.1068 2.8412 3.1798 3.2297 219 2.5597 3.3481 3.3487 3.3663 1.0855 2.9101 3.2766 3.3307

Within the data obtained by the different methods, we get that the approximate values obtained by Clark’s algorithm differ from the sample means at most by 57.6% when H=0.09 , by 13.08% when H=0.01 , by 2.85% when H=0.0013 , and by 1.06% when H=0.0001 . Thus, when H0.0013 , the values of the expected maximum, obtained by these completely different methods, are numerically identical. This indicates that the sample mean is approximately equal to Emaxi=1,NBH(i/N) .

Bounds for the approximation error

In this section, we find bounds for the error of approximation (2). As noted before, Emaxt[0,1]BH(t) (4Hπeln2)1/2 . It is expected that obtained sample means of the maximum functional (6) also satisfies this constraint. In Fig. 3, the sample means and the values of (4Hπeln2)1/2 are presented. As one can see, the inequality Emaxi=1,NBH(i/N)(4Hπeln2)1/2 is false for small values of H.

Sample means of the maximal functional

There are two possible explanations of this fact: either there is a significant error in calculations, or the approximation error ΔN grows rapidly as H0 . Let us verify these two explanations.

From [1, Theorem 4.2] we get that the expectation of the maximal functional (6) grows as H0 and has the limit limH0Emaxi=1,NBH(i/N)=12E(maxi=1,Nξi)+, where ξ1,,ξN are i.i.d. r.v.s, ξ1N(0,1) , and x+:=max{0,x} .

Moreover, the rate of convergence in (8) is also obtained in : 012E(maxi=1,Nξi)+Emaxi=1,NBH(i/N)11N2H. The right-hand side of (9) does not exceed 0.1 when N=220 and H<0.0038 .

We apply two approaches to calculate 12E(maxi=1,Nξi)+ . The first one is Monte Carlo simulations.

The sample means of 12(maxi=1,Nξi)+ are presented in Table 2 for several sample sizes n.

As we see, with increasing sample size 20 times, the sample means differ at most by 0.33% for each N. Therefore, to ensure the accuracy of calculations, it suffices to put n=1000 . Under such conditions, technical resources allow us to calculate the sample means for larger values of N. In Table 3, the values of the sample means are presented for N={220,221,222,223,224,225} .

Values of limit (8)

 Sample means of 12(maxi=1,N‾ξi)+ N∫0.51erf(−1)(2z−1)zN−1dz N\n 1000 5000 10000 15000 20000 28 1.9908 1.9908 1.9957 1.9965 1.9961 1.9989 29 2.1462 2.1506 2.1520 2.1526 2.1525 2.1524 210 2.3071 2.3033 2.3006 2.3004 2.2994 2.2969 211 2.4409 2.4362 2.4360 2.4371 2.4351 2.4337 212 2.5712 2.5657 2.5648 2.5635 2.5643 2.5640 213 2.6824 2.6847 2.6877 2.6874 2.6867 2.6887 214 2.8150 2.8066 2.8065 2.8060 2.8078 2.8082 215 2.9190 2.9259 2.9235 2.9248 2.9244 2.9232 216 3.0301 3.0372 3.0353 3.0340 3.0348 3.0343 217 3.1387 3.1372 3.1424 3.1418 3.1414 3.1417 218 3.2394 3.2456 3.2469 3.2460 3.2461 3.2458 219 3.3402 3.3442 3.3450 3.3458 3.3460 3.3469

Values of limit (8) for N220

 N 220 221 222 223 224 225 Sample means of 12(maxi=1,N‾ξi)+ 3.4516 3.536 3.627 3.724 3.816 4.073 N∫0.51erf(−1)(2z−1)zN−1dz 3.4452 3.541 3.634 3.726 3.815 3.902

Instead of generating random samples, we may calculate the value of 12E(maxi=1,Nξi)+ as an integral.

Let ξ1,,ξN be i.i.d. r.v.s, ξ1N(0,1) . Then 12E(maxi=1,Nξi)+=N21/21Φ(1)(z)zN1dz=N1/21erf(1)(2z1)zN1dz, where Φ(1) is the inverse function of Φ(x)=xey2/22πdy , xR , and erf(1) is the inverse function of the error function erf .

The proposition follows straightforwardly by quantile transformation.  □

We immediately get the following corollary.

For any H(0,1) and N1 , we have Emaxi=1,NBH(i/N)N1/21erf(1)(2z1)zN1dz.

The integrand in (11) is not an elementary function, but its values are tabulated, and there exist methods for its numerical computing. For the present paper, the integral N0.51erf(1)(2z1)zN1dz is calculated numerically, and the corresponding values are presented in Tables 2 and 3. By maintaining the accuracy of calculations, the maximum possible value of N is 231, and the value of the integral reaches 4.390.

The values of 12E(maxi=1,Nξi)+ , obtained by the two methods, differ at most by 0.44 % when N224 . When N=220 , the absolute error of numerical computing of (10) is less than 1.3×105 . Thereafter, for N=220 , inequality (11) becomes Emaxi=1,NBH(i/N)3.4452,H(0,1).

Let us return to the lower bound for Emaxi=1,NBH(i/N) . By Sudakov’s inequality  we have Emaxi=1,NBH(i/N)ln(N+1)N2H2πln2.

Moreover, the maximum of the right-hand side of (13) equals (4Hπeln2)1/2 and is reached when N=[e1/2H] . The values of the lower bound are presented in Table 4.

Lower bounds

 H 0.5000 0.0900 0.0100 0.0013 0.0001 (2Hπeln2)−1 0.5811 1.3696 4.1089 11.396 41.089 e1/2H 2.7183 258.67 5.18 × 1021 1.1 × 10167 2.97 × 102171 N (ln(N+1)N2H2πln2)1/2 28 0.0705 0.6853 1.0679 1.1207 1.1282 29 0.0529 0.6828 1.1246 1.1873 1.1963 210 0.0394 0.6761 1.1772 1.2503 1.2608 211 0.0292 0.6662 1.2260 1.3101 1.3222 212 0.0216 0.6537 1.2717 1.3671 1.3808 213 0.0159 0.6393 1.3145 1.4217 1.4371 214 0.0117 0.6233 1.3547 1.4740 1.4913 215 0.0085 0.6061 1.3925 1.5244 1.5435 216 0.0062 0.5881 1.4283 1.5729 1.5940 217 0.0045 0.5696 1.4620 1.6199 1.6429 218 0.0033 0.5507 1.4940 1.6653 1.6905 219 0.0024 0.5315 1.5244 1.7094 1.7367

Combining Tables 1, 2, and 4, we get that all obtained sample means for Emaxi=1,NBH(i/N) satisfy the constraint (ln(N+1)N2H2πln2)1/2Emaxi=1,NBH(i/N)N1/21erf(1)(2z1)zN1dz. Therefore, even with small values of the parameter H, the simulation does not lead to contradiction.

Now let us find a lower bound for the approximation error ΔN . We prove the following proposition.

Let ΔN be defined by (3). Then, for any H(0,1) and N1 , we have ΔN12Hπeln2N1/21erf(1)(2z1)zN1dz.

The statement follows from inequalities (1) and (11).  □

From this it follows that, for a fixed N, the approximation error ΔN+ as H0 . We also have the following evident corollaries.

Let N=220 . Then ΔN0.2055H3.4452,H(0,1).

The statement follows from inequalities (13) and (14).  □

Let N=220 . Then for the relative error, we have δH:=ΔNEmaxt[0,1]BH(t)116.765H,H(0,1).

The statement follows from inequalities (1) and (12).  □

When N=220 , from inequalities (15) and (16) we get the following conclusions:

if H<0.00022 , then the relative error δH75% , and ΔN>10.34 10.34$]]>; if H<0.00089 , then the relative error δH50% , and ΔN>3.45 3.45$]]>;

if H<0.0020 , then the relative error δH25% , and ΔN>1.15 1.15$]]>; if H<0.0028 , then the relative error δH10% , and ΔN>0.38 0.38$]]>;

if H<0.0032 , then the relative error δH5% , and ΔN>0.18 0.18$]]>; if H<0.0035 , then the relative error δH1% , and ΔN>0.03 0.03$]]>.

Thus, we conclude that the estimation of Emaxt[0,1]BH(t) by Monte Carlo simulations leads to significant errors for small values of the parameter H.

Acknowledgments

The author is grateful to prof. Yu. Mishura for numerous interesting discussions and active support.

References Borovkov, K., Mishura, Y., Novikov, A., Zhitlukhin, M.: Bounds for expected maxima of Gaussian processes and their discrete approximations. Stoch. Int. J. Probab. Stoch. Process. (2015). doi:10.1080/17442508.2015.1126282 Clark, C.E.: The greatest of a finite set of random variables. Oper. Res. 9(2), 145162 (1961). MR0125604 Coeurjolly, J.-F.: Simulation and identification of the fractional Brownian motion: A bibliographical and comparative study. J. Stat. Softw. 5(7), 153 (2000) Delorme, M., Wiese, K.J.: Maximum of a fractional Brownian motion: Analytic results from perturbation theory. Phys. Rev. Lett. 115(21), 210601 (2015) Molchan, G.M.: Maximum of a fractional Brownian motion: Probabilities of small values. Commun. Math. Phys. 205(1), 97111 (1999). MR1706900. doi:10.1007/s002200050669 Sinai, Y.G.: Distribution of the maximum of a fractional Brownian motion. Russ. Math. Surv. 52(2), 359378 (1997). MR1480141. doi:10.1070/RM1997v052n02ABEH001781 Sudakov, V.N.: Geometric Problems in the Theory of Infinite-Dimensional Probability Distributions vol. 141. Am. Math. Soc. (1979). MR0530375 Talagrand, M.: Lower classes for fractional Brownian motion. J. Theor. Probab. 9(1), 191213 (1996). MR1371076. doi:10.1007/BF02213740 Wood, A.T., Chan, G.: Simulation of stationary Gaussian processes in [0, 1]d. J. Comput. Graph. Stat. 3(4), 409432 (1994). MR1323050. doi:10.2307/1390903