Sampling first-passage times of fractional Brownian motion using adaptive bisections

We present an algorithm to efficiently sample first-passage times for fractional Brownian motion. To increase the resolution, an initial coarse lattice is successively refined close to the target, by adding exactly sampled midpoints, where the probability that they reach the target is non-negligible. Compared to a path of $N$ equally spaced points, the algorithm achieves the same numerical accuracy $N_{\mathrm{eff}}$, while sampling only a small fraction of all points. Though this induces a statistical error, the latter is bounded for each bridge, allowing us to bound the total error rate by a number of our choice, say $P_{\mathrm{error}}^{\mathrm{tot}}={10}^{-6}$. This leads to significant improvements in both memory and speed. For $H=0.33$ and $N_{\mathrm{eff}}=2^{32}$, we need $5000$ times less CPU time and $10000$ times less memory than the classical Davies-Harte algorithm. The gain grows for $H=0.25$ and $N_{\mathrm{eff}}=2^{42}$ to $3\times {10}^5$ for CPU and ${10}^6$ for memory. We estimate our algorithmic complexity as ${\mathcal{C}}^{\mathrm{ABSec}}\left(N_{\mathrm{eff}}\right)=O\left[{\left(ln{N}_{\mathrm{eff}}\right)}^3\right]$, to be compared to Davies-Harte, which has complexity ${\mathcal{C}}^{\mathrm{DH}}\left(N\right)=O\left(NlnN\right)$. Decreasing $P_{\mathrm{error}}^{\mathrm{tot}}$ results in a small increase in complexity, proportional to $ln(1/P_{\mathrm{error}}^{\mathrm{tot}})$. Our current implementation is limited to the values of $N_{\mathrm{eff}}$ given above, due to a loss of floating-point precision. Our algorithm can be adapted to other extreme events and arbitrary Gaussian processes. It enables one to numerically validate theoretical predictions that were hitherto inaccessible.

Figure 1

The continuous stochastic path (gray rough line) crosses the barrier (blue horizontal line) for the first time at ${\tau }^{\infty }$ (black leftmost square mark). The discretization with $N$ points (red line passing through rightmost square) overestimates this time as ${\tau }^N$ (red rightmost square mark). The numerical estimate is improved to ${\tau }^{4N}$ (green middle square mark) when refining the discretization (green line passing through middle square mark). This systematic error worsens for diminishing values of Hurst parameter $H$.

Figure 2

Illustration of the adaptive bisection routine. The grid $\mathcal{T}$ (bottom) contains points in time; here details are shown of the initial bridge $t_{\mathrm{l}}=i{2}^{-g},t_{\mathrm{r}}=(i+1)2^{-g}$ (labeled bullets) and successively introduced midpoints (bullets on the time axis). The path $\mathcal{X}$ (above) samples values at times (dashed lines) which approximate the path by linear interpolations (gray and black thick lines). The threshold $m$ (red uppermost horizontal line) is crossed by the path and bisections are generated for every bridge the end points of which lie in the critical strip corresponding to its level (blue vertical lines underneath). The horizontal arrows on top of the path indicate the bridges in between the grid points. The mapping from bridges to the binary tree (top) is indicated with dotted lines. The top node (1) corresponds to the widest bridge $(i{2}^{-g},(i+1)2^{-g})$, and children correspond to subintervals generated by the midpoint. The bridges are explored in order as given by numbers above nodes and chosen by the bridge-selection routine (see text for details). Bridges that are critical (blue filled nodes) are bisected, and their children are checked from left to right, until a first-passage event has been identified at maximum bisection level $L$ (red filled node 7). This event terminates the algorithm. In contrast to node 1, which belongs to the initial grid ${\mathcal{T}}^{\left(0\right)}$, nodes 2 to 7 stem from adaptive bisections and contribute to the total count of bisections $M$. The maximum number of nodes which could theoretically be spawned off this particular subinterval is $2^{L-g}-1$.

Figure 3

Error rate from the phone book test for various values of ${\varepsilon }^{\prime }$ for $H=0.33$. The error rate is almost identical when changing the initial grid size from $2^8$ (red square marks) to $2^4$ (blue triangle marks) at the same maximum bisection level $L=20$. When lowering the maximum bisection level to $L=16$, the error rate improves. The relation between error rate and error tolerance decreases approximatively linearly over several orders of magnitude (compare with solid gray line). The total error rate is approximately $10{\varepsilon }^{\prime }$ for $L=20$ (solid gray line) and about $3{\varepsilon }^{\prime }$ for $L=16$ (dashed gray line). Note that the prefactor is much smaller than the number of points, which can be read off from Fig. 5. Error rates were averaged over ${10}^5$ to ${10}^6$ iterations.

Figure 4

Average number of new midpoints generated at bridge level $\ell$, for various values of $H$ (solid, dashed, dash-dotted, and dotted lines) as a function of $\ell H$. For equal values of $\ell H$, a lower Hurst parameter implies a larger number of average bisections. These numbers are virtually independent of the initial grid size, as is shown for ${\mathrm{\Lambda }}^4$ (circle marks) and ${\mathrm{\Lambda }}^8$ (triangle marks).

Figure 5

Average number of bisections $M$ as a function of the maximum bisection level $L$ (i.e., $N_{\mathrm{eff}}=2^L$) for different values of $H$ (diamond, circle, upright, and upside down triangle marks). Inset shows $M$ vs $LH$. As long as $H\ge 0.33$ growth is asymptotically approximately linear in $L$, corroborating $M\sim ln\left(N_{\mathrm{eff}}\right)$. For smaller values of $H$, either the linear regime is not yet reached, or the growth is stronger (5000 iterations with initial grid ${\mathrm{\Lambda }}^8$ and error tolerance ${\varepsilon }^{\prime }={10}^{-9}$). For $H=0.5$, extrapolation was used.

Figure 6

Average user time required to find the first-passage time in a grid of effective discretization precision $2^{-LH}$. The dashed lines indicate user time for Davies-Harte method, and solid lines are for the adaptive bisection method. The three different colors indicate $H=0.33,0.5,0.67$ (bottom, center, top pairs of lines). Simulations were run ${10}^4$ times for ${\varepsilon }^{\prime }={10}^{-9}$ and for two different initial subgrid sizes (${\mathrm{\Lambda }}^4$, circle marks; ${\mathrm{\Lambda }}^8$, square marks). For $H=0.33$ (top solid blue lines), the effective system sizes range from $L=4$ to 32 for ${\mathrm{\Lambda }}^4$, and $L=12$ to 28 for ${\mathrm{\Lambda }}^8$. For $H=0.5$ (center solid green lines), $L$ ranges from 4 to 22 for ${\mathrm{\Lambda }}^4$ and from $L=12$ to 22 for ${\mathrm{\Lambda }}^8$. For $H=0.67$ (bottom solid red lines), $L$ ranges from 4 to 16 for ${\mathrm{\Lambda }}^4$ and 12 to 16 for ${\mathrm{\Lambda }}^8$.

Figure 7

${(t_{\mathrm{user}}/\mathrm{iteration})}^{1/3}$ plotted vs effective discretization $N^H$ for various values of $H$ (blue circle marks $H=0.33$, green square marks $H=0.5$, red diamond marks $H=0.67$; see Fig. 6). They corroborate the estimate of ${\mathcal{C}}^{\mathrm{ABSec}}\sim {\left(ln{N}_{\mathrm{eff}}\right)}^3$. Straight lines indicate fits of the form $aln\left(N\right)+b$ vs $t_{\mathrm{user}}^{1/3}$ implying a scaling of $t_{\mathrm{user}}\sim a^3\left[ln{\left(N_{\mathrm{eff}}\right)}^3\right]+O\left[{\left(ln{N}_{\mathrm{eff}}\right)}^2\right]$. This is in agreement with the complexity estimate in Eq. (35). The inset shows the ratio between data points and the fit.

Figure 8

User time for ABSec (solid lines) compared to DH (dashed line) for two different initial grid sizes and two different values of error tolerance ${\varepsilon }^{\prime }$. For a 100 times higher error tolerance (top semitransparent pair of lines), user time increases by up to $60\%$.

Figure 9

Memory usage for DH (dashed line) and ABSec (solid line) for two different initial subgrid sizes. DH scales linearly in $N$, while ABSec grows only slowly (see text for estimate). For system of size $N_{\mathrm{eff}}=2^{28}$, ABSec needs only ${10}^{-2}$ to ${10}^{-3}$ of the memory for DH. For larger systems or smaller $H$, the advantage of ABSec is even bigger. Measurements were taken after ${10}^4$ iterations.

Figure 10

Ratio between the sampled variance and no-neighbor-estimate of variance [see Eq. (17)] of an inserted midpoint $X_m$ vs the level of the bisected bridge. For $H=0.5$ (green diamond marks), the ratio equals 1, as BM is Markovian. For $H\ne 0.5$ (red square marks $H=0.67$, blue circle lines $H=0.33$), the variance fluctuates, as shown by the error bars for one standard deviation. Numerical errors due to a loss of floating-point precision become relevant around $L_{\max }\simeq 11/H$. ABSec was used with an initial grid ${\mathrm{\Lambda }}^8$ and ${\varepsilon }^{\prime }={10}^{-9}$.