A fractional Brownian motion {BH(t),t≥0} is a centered Gaussian stochastic process with covariance function E[BH(t)BH(u)]=12(t2H+u2H−|t−u|2H),t,u≥0, 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),t≥0}=d{a−HBH(at),t≥0}, 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(T−1t),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),t≥0} . The distribution of maxt∈[0,1]W(t) is known. Namely, P(maxt∈[0,1]W(t)≤x)=2π∫0xe−y2/2dy,x≥0, 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 [5] has found an asymptotic behavior of small-ball probabilities for the maximum of the fractional Brownian motion. Talagrand [8] 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 [6] and recently by Delorme and Weise [4].

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. [1] we know the following bounds: (1) 12Hπeln2≤Emaxt∈[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, (2) Emaxt∈[0,1]BH(t)≈Emaxi=1,N‾BH(i/N). The authors of [1] obtain an upper bound for the error ΔN of approximation (2). Namely, for N≥21/H , (3) 0≤ΔN:=Emaxt∈[0,1]BH(t)−Emaxi=1,N‾BH(i/N) (4) ≤2lnNNH(1+4NH+0.0074(lnN)3/2).

The implementation of approximation (2) has technical limitations. Due to modern computer capabilities, we assume that N≤220≈106 . Under such conditions, inequality (4) is true when H≥0.05 , and ΔN<11.18 .

In this article, we make Monte Carlo simulations and estimate Emaxi=1,N‾BH(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 H→0 . For the rate of this convergence, when N=220 , we prove the inequality ΔN>c1H−1/2−c2 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 .

In this section, we present results of approximate computing Emaxi=1,N‾BH(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 {10−4(1+4i),i=0,24‾}∪{10−2i,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: (6) maxi=1,N‾BH(i/N), (7) 1N ∑i=1NBH(i/N).

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,N‾BH(i/N)≥(4Hπeln2)−1/2 is false for small values of H.

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

From [1, Theorem 4.2] we get that the expectation of the maximal functional (6) grows as H→0 and has the limit (8) limH→0Emaxi=1,N‾BH(i/N)=12E(maxi=1,N‾ξi)+, where ξ1,…,ξN are i.i.d. r.v.s, ξ1∼N(0,1) , and x+:=max{0,x} .

Moreover, the rate of convergence in (8) is also obtained in [1]: (9) 0≤12E(maxi=1,N‾ξi)+−Emaxi=1,N‾BH(i/N)≤1−1N2H. 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} .

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

We immediately get the following corollary.

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 N∫0.51erf(−1)(2z−1)zN−1dz 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 2³¹, 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 N≤224 . When N=220 , the absolute error of numerical computing of (10) is less than 1.3×10−5 . Thereafter, for N=220 , inequality (11) becomes (12) Emaxi=1,N‾BH(i/N)≤3.4452,H∈(0,1).

Let us return to the lower bound for Emaxi=1,N‾BH(i/N) . By Sudakov’s inequality [1, 7] we have (13) Emaxi=1,N‾BH(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.

Combining Tables 1, 2, and 4, we get that all obtained sample means for Emaxi=1,N‾BH(i/N) satisfy the constraint (ln(N+1)N2H2πln2)1/2≤Emaxi=1,N‾BH(i/N)≤N∫1/21erf(−1)(2z−1)zN−1dz. 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.

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

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

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.