Adaptive Quantile Computation for Brownian Bridge in Change-Point Analysis
Jürgen Franke, Mario Hefter, André Herzwurm, Klaus Ritter, Stefanie Schwaar
AADAPTIVE QUANTILE COMPUTATION FOR BROWNIANBRIDGE IN CHANGE-POINT ANALYSIS
JÜRGEN FRANKE, MARIO HEFTER, ANDRÉ HERZWURM, KLAUS RITTER,AND STEFANIE SCHWAAR
Abstract.
As an example for the fast calculation of distributional parametersof Gaussian processes, we propose a new Monte Carlo algorithm for the com-putation of quantiles of the supremum norm of weighted Brownian bridges. Asit is known, the corresponding distributions arise asymptotically for weightedCUSUM statistics for change-point detection. The new algorithm employs anadaptive (sequential) time discretization for the trajectories of the Brownianbridge. A simulation study shows that the new algorithm by far outperformsthe standard approach, which employs a uniform time discretization. Introduction
In statistical inference, asymptotics frequently leads to the distribution of nonlin-ear functionals of Gaussian processes. E.g., the construction of uniform asymptoticconfidence bands for a regression function based on kernel estimates requires thestudy of the supremum of the absolute values of a certain Gaussian process, cf.Härdle (1990, Sec. 4.3). To mention another example, for testing the equality ofmean functions in functional data analysis, a test statistic is used which underthe hypothesis converges in distribution to an integral of the square of a Gaussianprocess, cf. Horváth and Kokoszka (2012, Sec. 5.1).In some cases like the first example above, cf. Bickel and Rosenblatt (1973), thedistribution of the nonlinear functional may be derived in a form which, in par-ticular, allows for the calculation of quantiles for tests and confidence assessments.In many other cases, however, distributional characteristics have to be calculatednumerically by Monte Carlo simulation. Even in rather simple cases where theGaussian process is just a Wiener process or a Brownian bridge, using the standardapproximation of a continuous-time Gaussian process by a corresponding discrete-time process on an equidistant grid may result in a severe computational load if adecent quality of approximation is required. We shall discuss this below in moredetail for a specific case.In this paper, we show that this problem can be overcome by using fast adaptiveapproximation methods for the strong (or pathwise) approximation of nonlinearfunctionals of Gaussian processes. Adaptive algorithms employ sequential strategiesto construct the underlying discretization, which may therefore be adjusted to keyfeatures of the individual trajectories.For illustrating our approach, we focus on weighted CUSUM tests for change-points, which leads to the distribution of the supremum of a weighted reflectingBrownian bridge, i.e., the supremum norm of a weighted Brownian bridge. Westress that the basic idea can be transferred to other situations where, e.g., quantilesof nonlinear functionals of Gaussian processes have to be calculated by Monte Carlosimulation.
Date : November 12, 2020.
Key words and phrases. change-point problem, weighted CUSUM statistic, weighted Brown-ian bridge, sup-norm quantiles, Monte Carlo algorithm, adaptive discretization. a r X i v : . [ s t a t . C O ] D ec J. FRANKE, M. HEFTER, A. HERZWURM, K. RITTER, AND S. SCHWAAR
Change-point tests are of interest in many areas of applications, e.g., in produc-tion monitoring, see Page (1957), on-line-monitoring of intensive-care patients, seeFried and Imhoff (2004), or global warming studies, see Gallagher, Lund, and Rob-bins (2013), to name just a few. The first publications about testing for a changein data go back to the 1950s, see, e.g., Page (1957), who has considered testing fora change in the mean and has used weighted cumulated sums of sample residuals,so called weighted CUSUM statistics. The corresponding weight function is givenby w ( t ) = ( t · (1 − t )) − / for 0 < t <
1, so that the weighted cumulated sums are variance stable. The distri-bution of those statistics is determined asymptotically in a variant of the Darling-Erdös Theorem, see Theorem 1, which immediately yields asymptotic quantiles.To cope with performance problems regarding size and power, the standardweights of CUSUM statistics have been modified which results asymptotically inthe distribution of the supremum of a weighted reflecting Brownian bridge. Thosestatistics have no size problems and better power against changes of the mean closerto the boundaries of the observation interval. Simulation studies using differentweight functions and analyzing the power of the corresponding CUSUM-type testsfor different positions of the change (early, middle and late) show the importance ofthe weight function, see Csörgő and Horváth (1997), Kirch and Tadjuidje Kamgaing(2016), and Schwaar (2016). An overview to such general CUSUM-type tests isgiven in Aue and Horváth (2013).In the present paper we consider weight functions of the form w η,γ ( t ) = 1 ] η, − η [ ( t ) · ( t · (1 − t )) − γ for 0 < t < ≤ η < / ≤ γ ≤ /
2, and the correspondingconvergence result for ( η, γ ) = (0 , /
2) is formulated in Theorem 2. Except forthe extremal cases γ = 0 and γ = 1 /
2, there is no known method for analyticallycalculating quantiles or other characteristics of the limit distributions for these ormore general weight functions. Hence, we have to use Monte Carlo simulationwhere a crucial part consists in generating paths of a Brownian bridge. This is inparticular computationally very expensive if we are interested in calculating, e.g.,extreme quantiles beyond the common levels 0 .
05 or 0 .
01 with a high precision upto 10 − . Such extreme levels of confidence are common in many applications inindustry or medicine where a high degree of reliability is required. Also, in viewof the Bonferroni inequality, low p -values of tests are of interest in multiple testingsituations including many hypotheses, see, e.g., Hochberg (1988).In this paper, we propose an adaptive algorithm which reduces the computationtime for calculating asymptotic quantiles of CUSUM test statistics with weightfunction w η,γ for γ
6∈ { , / } , where the reduction turns out to be dramatic inthe more challenging situations. Consider, for instance, the task to compute the0 . − for the parameters η = 0 and γ = 0 .
45. Ona common up-to-date processor the standard algorithm with an equidistant gridrequires more than two hours of computation time, while the adaptive algorithmachieves the same goal within 12 seconds. The reason for this is, roughly speaking,that the adaptive algorithm allows to sample almost exactly from the correct limitdistribution at a reasonable computational cost.This paper is organized as follows. In Section 2 we give a brief sketch of thechange-point application that is used as a demonstrator for our approach. InSections 3 and 4 we consider the strong approximation of the supremum of the
DAPTIVE QUANTILE COMPUTATION FOR BROWNIAN BRIDGE 3 unweighted Brownian motion and of weighted reflecting Brownian bridges, respec-tively. For the former problem, Calvin, Hefter, and Herzwurm (2017) have con-structed an adaptive algorithm that strikingly outperforms all non-adaptive al-gorithms, see Theorems 7 and 9. See also Calvin (1997, 2001, 2004) for relatedconvergence results. No such result is known for weighted reflecting Brownianbridges, but still we construct a modification of the adaptive algorithm for theseprocesses in Section 4.1. Numerical experiments reveal again the vast superiorityof the adaptive algorithm over, at least, the standard algorithm that is based onan equidistant grid, see Section 4.2. In Section 5 we study quantile computationfor the supremum of a weighted reflecting Brownian bridge. We present a newalgorithm, with the adaptive algorithm from Section 4.1 as the key ingredient, thatyields a quantile up to a user-specified error tolerance, see Section 5.1. Numericalexperiments, which show the superiority of the new algorithm over the standardapproach, are presented in Section 5.2.2.
The Statistical Problem: Change-Point Test
We are interested in detecting a structural change, specifically at most one change(AMOC), in a time series model, and for illustration we consider the following mostsimple model with a possible mean change. The model, which is one of the earliestchange-point models analyzed, is given by ξ i = ( ε i , if i ≤ m,d + ε i , if i > m, for i = 1 , . . . , n , where n, m ∈ N with n ≥ ≤ m = m ( n ) ≤ n, and where d ∈ R with d = 0, see Page (1957). The residuals ε i are assumed to beiid, each centered with finite second moment σ >
0, which is assumed to be knownfor simplicity. If m < n , then a structural change is present, and m is called thechange-point. A test is constructed for H : m = n, H : m < n. Based on the quasi likelihood ratio test, the weighted CUSUM statistic T n ( w ) = max ≤ k Theorem 1 (Darling-Erdös Theorem) . Let w ( t ) = ( t · (1 − t )) − / for < t < , and assume that E( | ε i | δ ) < ∞ for some δ > . Under H we have lim n →∞ P ( { T n ( w ) ≤ c n ( α ) } ) = 1 − α J. FRANKE, M. HEFTER, A. HERZWURM, K. RITTER, AND S. SCHWAAR for < α < and c n ( α ) = σa (log n ) (cid:0) − log (cid:0) − log(1 − α ) (cid:1) + b (log n ) (cid:1) , with a ( x ) = p x and b ( x ) = 2 log x + log log x − log π. Theorem 1 immediately yields an asymptotic level α test, see Remark 4. How-ever, for small sample sizes n the convergence in the Darling-Erdös Theorem oftenleads to level distortion, see Kirch (2006). To overcome this problem modificationsof the weight function w are considered, see, e.g., Csörgő and Horváth (1988). Inthis paper we study weight functions w η,γ of the form(1) w η,γ ( t ) = 1 ] η, − η [ ( t ) · ( t · (1 − t )) − γ for 0 < t < 1, where 0 ≤ η < / , ≤ γ ≤ / . Observe that Theorem 1 deals with the case ( η, γ ) = (0 , / X = ( X ( t )) t ∈ ]0 , we put S ( X ) = sup Let w η,γ be given by (1) with ( η, γ ) = (0 , / . Under H we have σ · T n ( w η,γ ) d −→ S ( w η,γ | B | ) . Basically, Theorem 2 yields an asymptotic level α test. However, for applica-tion the quantiles of the supremum S ( w η,γ | B | ) of the weighted reflecting Brownianbridge w η,γ | B | are needed, i.e., they have to be known analytically or to be eas-ily computed numerically, see Remark 4. In Schwaar (2020) besides data drivenweighted change-point estimators, data driven weighted change-point tests are con-sidered, where the parameter γ is replaced by an estimator and η = 0. Especiallythen knowledge about the quantiles for the weighted test statistic with non-extremevalues is required.For completeness we add that S ( w , / | B | ) = ∞ with probability one, whichfollows from the law of the iterated logarithm, and hence Theorem 2 implies T n ( w , / ) p −→ ∞ . Remark 3. Consider the extremal cases γ = 0 and γ = 1 / 2. For ( η, γ ) = (0 , γ = 0 and 0 ≤ η < / γ = 1 / < η < / S ( w η, / | B | ) and related quantities, see also Yor (1984). The Mellin transform of acorresponding hitting time for a standard Brownian motion has been determined,and values of the c.d.f. have been obtained via numerical inversion. DAPTIVE QUANTILE COMPUTATION FOR BROWNIAN BRIDGE 5 To the best knowledge of the authors, no (semi-)analytic way to determine thequantiles is known beyond the extremal cases, i.e., for 0 < γ < / ≤ η < / Remark 4. Let c η,γ,n > 0. Consider the test that rejects H if and only if | T k,n |√ n > c η,γ,n · (cid:18) kn · (cid:18) − kn (cid:19)(cid:19) γ for some k ∈ N with η · n < k < (1 − η ) · n. To obtain an asymptotic level α test we proceed as follows. For 0 ≤ γ < / c η,γ = c η,γ,n is determined by P ( { S ( w η,γ | B | ) ≤ c η,γ /σ } ) = 1 − α, see Theorem 2, which may be solved semi-analytically if γ = 0 and numerically,using the algorithm presented in Section 5, if 0 < γ < / 2. For γ = 1 / c η,γ,n . Example 5. We illustrate the role of the parameter γ in the case η = 0. Let σ = 1and α = 0 . 05. In Figure 1 the threshold function f γ,n ( t ) = c ,γ,n · ( t · (1 − t )) γ is presented for γ = 0 , . , . 45, and for γ = 0 . n = 10 , . The criticalvalues are given numerically as follows. Remark 3 yields c , . = 1 . Q ad ε ( w, q ) according to Section 5 with ε = 10 − yields c , . =1 . 99 and c , . = 2 . 91. Furthermore, c , . , ≈ . 241 and c , . , ≈ . 353 dueto Theorem 1. 0 0 . . . t f γ , n ( t ) γ = 0 γ = 0 . γ = 0 . γ = 0 . n = 10 γ = 0 . n = 10 Figure 1. Threshold function f γ,n from Example 5 for η = 0, σ = 1, and α = 0 . γ close to 1 / γ = 0 has a higher power than the one with γ close to1 / J. FRANKE, M. HEFTER, A. HERZWURM, K. RITTER, AND S. SCHWAAR Remark 6. The construction of change-point tests in more complicated time seriesmodels, e.g., with a serial dependence, may also be based on CUSUM statisticsand quantiles of S ( X ) for suitable processes X . See Aue and Horváth (2013) formodels that lead to X = w | B | with weight functions w , or to X = P ‘j =1 B j withindependent standard Brownian bridges B j . More generally, see Aue, Rice, andSönmez (2018) for functional time series models that lead to X = P ∞ j =1 λ j B j withnon-negative scalars λ j . Our approach in Section 5 can easily be adapted to theseclasses of processes.3. Approximation of the Supremum of a Brownian Motion In this section we discuss the strong (or pathwise) approximation of the supre-mum S ( W ) of a standard Brownian motion W = ( W ( t )) t ∈ [0 , on the unit interval.We consider algorithms A that evaluate W at a finite number of points t k ∈ ]0 , S ( W ) by the discrete maximum of W at these points. The error e ( A ) of any such measurable algorithm A is defined by e ( A ) = E ( | S ( W ) − A ( W ) | ) . We recall known results that demonstrate that suitable adaptive algorithms, i.e.,algorithms that sequentially evaluate any trajectory of W , are far superior to non-adaptive algorithms, i.e., algorithms that are based on a fixed, a priori given dis-cretization of ]0 , n evaluationsof W . The following result is due to Ritter (1990). Theorem 7. There exist constants c , c > with the following properties for every n ∈ N . The algorithm A eq n given by A eq n ( W ) = max k =1 ,...,n W ( k/n ) satisfies (2) e ( A eq n ) ≤ c · n − / . For every choice of t , . . . , t n ∈ ]0 , the algorithm A non n given by A non n ( W ) = max k =1 ,...,n W ( t k ) satisfies (3) e ( A non n ) ≥ c · n − / . For the purpose of this paper the crucial part of Theorem 7 is the lower bound (3),which says that non-adaptive algorithms achieve at most the order of convergence1 / 2. This lower bound is sharp, up to constants, as we have a matching upperbound (2), which is already achieved by an equidistant discretization. Remark 8. We conjecture that the lower bound (3) is also valid for algorithms A non n of the form A non n ( W ) = ζ ( W ( t ) , . . . , W ( t n )), where ζ : R n → R is any measurablemapping. The conjecture is true at least if t k = k/n , see Hefter and Herzwurm(2017, Prop. 2.1 and proof of Thm. 3.3), or if the L -error of A non n is considered,instead of the L -error e ( A non n ), see Calvin (2004, beginning of proof of Thm. 2.1).The asymptotic distribution of the error S ( W ) − A eq n ( W ) has been derived inAsmussen, Glynn, and Pitman (1995, Thm. 1).An adaptive algorithm A ad n that uses n evaluations of the Brownian motion W is formally defined by a point t ∈ ]0 , 1] and Borel-measurable mappings χ k : R k − → ]0 , DAPTIVE QUANTILE COMPUTATION FOR BROWNIAN BRIDGE 7 for k = 2 , . . . , n . Iteratively the algorithm computes y = W ( t ) and y k = W ( χ k ( y , . . . , y k − ))for k = 2 , . . . , n , and it yields the output A ad n ( W ) = max k =1 ,...,n y k . In this way the k -th evaluation site χ k ( y , . . . , y k − ) may depend on the previouslyobtained values y , . . . , y k − (and the corresponding evaluation sites). Of course,the non-adaptive algorithms, which are considered in Theorem 7, correspond to theparticular case of constant mappings χ k .The following result is due to Calvin, Hefter, and Herzwurm (2017). Theorem 9. There exists a sequence of adaptive algorithms A ad n with the followingproperty. For every ρ > there exists a constant c > such that for every n ∈ N e ( A ad n ) ≤ c · n − ρ . We refer to Calvin, Hefter, and Herzwurm (2017) for the construction of thealgorithms A ad n that are considered in Theorem 9; see also Section 4 for basic ideas.According to Theorem 9 suitable adaptive algorithms achieve, roughly speaking,the polynomial order of convergence ∞ . Combining Theorems 7 and 9 we seethat adaptive algorithms strikingly outperform all non-adaptive algorithms for thestrong approximation of the supremum S ( W ) of a Brownian motion.Of course, Theorems 7 and 9 are irrelevant for the computation of quantiles of S ( W ), since the distribution of S ( W ) is known explicitly. The theorems stronglysuggest, however, that adaptive algorithms should be considered for quantile com-putation if (semi-)analytic methods are not available. The latter holds true for theprocesses w η,γ | B | with γ 6∈ { , / } , see Section 2.4. Approximation of the Supremum of a Weighted ReflectingBrownian Bridge For notational convenience we put w = w η,γ , where ( η, γ ) = (0 , / The Adaptive Algorithm. In this section we present an adaptive algorithmfor the strong approximation of the supremum S ( w | B | ) of the weighted reflectingBrownian bridge w | B | . This algorithm is a modification of the algorithm con-structed in Calvin, Hefter, and Herzwurm (2017), which achieves the error boundin Theorem 9 for the strong approximation of the supremum S ( W ) of a standardBrownian motion W .Both of these adaptive algorithms are greedy algorithms, and the basic idea inthe construction is as follows. After the k -th step the algorithm has computeda partition of the interval [0 , 1] into k subintervals together with the values ofthe underlying stochastic process at the boundary points of all these intervals. Ascore value is available for each subinterval, and the subinterval with the highestscore will be split at the midpoint. Ideally, the score value should be the condi-tional probability that the corresponding subinterval contains a global maximizerof w | B | . Reasonable substitutes for these conditional probabilities are needed inthe computation.In the sequel we present the algorithm for the strong approximation of S ( w | B | )in detail. Based on the weight function w we assign a weight v ( ‘, r ) to any interval J. FRANKE, M. HEFTER, A. HERZWURM, K. RITTER, AND S. SCHWAAR [ ‘, r ] ⊆ [0 , 1] with positive length in the following way. Let c = ‘ + r . For a weight function w that is positive and differentiable with a ‘small’ derivativeeverywhere on ]0 , 1[ it is reasonable to take v ( ‘, r ) = w ( c ). Since these conditions arenot met for w = w η,γ , except for the trivial case ( η, γ ) = 0, we proceed differently.To avoid a too small score value we take v ( ‘, r ) = ( , if [ ‘, r ] ⊆ [0 , η ] ∪ [1 − η, , ( c · (1 − c )) − γ , otherwise.In fact, observe that v ( ‘, r ) = w ( c ) if and only if c ≤ η < r or ‘ < − η ≤ c . If v ( ‘, r ) = w ( c ), then w ( c ) = 0 while v ( ‘, r ) may potentially be very large.Suppose that x = B ( ‘ ) and y = B ( r ) are known, while no values of B are knowninside of ] ‘, r [. It is reasonable to compare v ( ‘, r ) · B ( c ) with a certain threshold m ,e.g., with the largest value of w | B | known so far, under the conditional distributionof B ( c ). The latter is the normal distribution with mean ( x + y ) / r − ‘ ) / 4. More precisely, we define the score function ϕ : D × R × [0 , ∞ [ → [0 , ∞ [by ϕ ( ‘, r, x, y, m ) = E (cid:0) ( v ( ‘, r ) · Z − m ) + (cid:1) + E (cid:0) ( v ( ‘, r ) · Z + m ) − (cid:1) , where D = { ( ‘, r ) ∈ [0 , : ‘ < r } and Z ∼ N (( x + y ) / , ( r − ‘ ) / . Remark 10. The score function is easily computed as follows. Let Φ denote thedistribution function of Y ∼ N (0 , a ∈ R . Put ψ ( a ) = Φ ( a ) + a · Φ( a ) . Since E (( Y + a ) + ) = ψ ( a ), we obtain ϕ ( ‘, r, x, y, m )= v ( ‘, r ) · √ r − ‘ · (cid:18) ψ (cid:18) x + y − m/v ( ‘, r ) √ r − ‘ (cid:19) + ψ (cid:18) − x + y + 2 m/v ( ‘, r ) √ r − ‘ (cid:19)(cid:19) if v ( ‘, r ) > 0. Otherwise ϕ ( ‘, r, x, y, m ) = 0, since m ≥ ψ one may uselim a →∞ ψ ( a ) a = 1and lim a →−∞ ψ ( a ) a − · Φ ( a ) = 1 . We use this asymptotic behavior for | a | > 3, i.e., we replace ψ by e ψ ( a ) = a, if a > ,a − · Φ ( a ) , if a < − ,ψ ( a ) , otherwise , in the computation of ϕ ( ‘, r, x, y, m ). DAPTIVE QUANTILE COMPUTATION FOR BROWNIAN BRIDGE 9 We are ready to define the adaptive algorithm, which will be denoted by A ad n ( w, · ).The algorithm will sequentially evaluate any trajectory of B , and the relevant in-formation about the corresponding partition of [0 , 1] after k steps is represented bya set I k as follows. There are k elements in I k , which correspond to the k subinter-vals in this partition. More precisely, ( ‘, r, x, y, s ) ∈ I k represents a subinterval [ ‘, r ]with boundary values x = B ( ‘ ) and y = B ( r ) and with score value s . Furthermore, m k ≥ w | B | after k steps. In the first step weput m = 0 , I = { (0 , , , , ϕ (0 , , , , m )) } . In the k -th step with 2 ≤ k ≤ n we choose any I = ( ‘, r, x, y, s ) ∈ I k − with themaximal value of s among all elements of I k − , i.e., with the largest score value(this choice needs not to be unique), and we evaluate B at the midpoint of thecorresponding subinterval, i.e., we compute z = B ( c ) . The new discrete maximum of w | B | is given by m k = max ( m k − , w ( c ) · | z | ) , and the new partition is represented by I k = ( I k − \ { I } ) ∪ { I , I } with I = ( ‘, c, x, z, ϕ ( ‘, c, x, z, m k )) , I = ( c, r, z, y, ϕ ( c, r, z, y, m k )) . After n steps the adaptive algorithm returns the output A ad n ( w, B ) = m n . Remark 11. We discuss the computational cost to simulate A ad n ( w, B ). First ofall, we use a max-priority queue for storing the elements of the sets I , . . . , I n . Thisallows to choose ( ‘, r, x, y, s ) ∈ I k − with the maximal value of s at a cost of orderln( k ) in the k -th step, see Cormen, Leiserson, Rivest, and Stein (2003, Sec. 6.5).Thereafter B ( c ) may be simulated at a constant cost, independently of k , under theconditional distribution. By definition, no updates of the score values are computedfor the elements in I k ∩ I k − , although m k may be larger than m k − . It followsthat the total cost to simulate A ad n ( w, B ) is of the order n ln( n ).4.2. Numerical Experiments. We compare the adaptive algorithm A ad n ( w, · ) andthe non-adaptive algorithm A eq n ( w, · ) given by A eq n ( w, B ) = max k =1 ,...,n − w ( t k ) · | B ( t k ) | with t k = k/n , which has been used by, e.g., Eastwood and Eastwood (1998) andOrasch and Pouliot (2004), and for a similar problem by Akashi, Dette, and Liu(2018).Analogously to Section 3 we consider the error(4) e ( A n ( w, · )) = E ( | S ( w | B | ) − A n ( w, B ) | )for A n ( w, · ) = A ad n ( w, · ) and A n ( w, · ) = A eq n ( w, · ). We add that exact values orupper or lower bounds of the error are not available for any of these algorithms.Therefore we determine the error via a Monte Carlo simulation, where we replace S ( w | B | ) in (4) by A ad n ( w, B ) or A eq n ( w, B ), respectively, with n being 10 timeslarger than the largest value of n that is considered in the numerical experiment.Simultaneously, we determine the average run-time of both algorithms. For eachvalue of n we use 10 Monte Carlo replications. − − − − n t i m e i n s ec γ = 0 . 45 (eq) γ = 0 . 25 (eq) γ = 0 . 45 (ad) γ = 0 . 25 (ad)order n order n ln( n ) Figure 2. Average run-time vs. number n of discretization pointsfor the strong approximation of S ( w | B | ).All programs are written in C++, and the computations are performed on asingle Intel Xeon Gold 6126 processor. All results are presented together withasymptotic confidence intervals with confidence level 0 . 95. In the numerical exper-iments we consider η = 0 and γ = 0 . 25 or γ = 0 . 45, cf. Example 5.In Figure 2 we relate the average run-time to the number n of discretizationpoints. For n ≤ the computational overhead due to adaption is moderate, asthe non-adaptive algorithm is at most 10 times faster than the adaptive algorithmin this range. Furthermore, the average run-times for both algorithms are in linewith the worst case behavior, namely, order n ln( n ) for A ad n ( w, · ), see Remark 11,and order n for A eq n ( w, · ). Finally, we see that the average run-time does not dependon γ .The relation between the error and the average run-time, which is most im-portant, is presented in Figure 3. First of all, we observe a polynomial order ofconvergence ∞ for the adaptive algorithm, in contrast to a polynomial order ofconvergence of only about 1 / γ , deteriorates the speed of convergence forboth algorithms. For the same average run-time of 10 − seconds the adaptive al-gorithm achieves an error of less than 10 − for both values of γ , while the errorof the non-adaptive algorithm is about 10 − . For quantile computation we maytherefore sample almost from the correct distribution by means of the adaptivealgorithm with a reasonable average run-time, while this is impossible by means ofthe non-adaptive algorithm. DAPTIVE QUANTILE COMPUTATION FOR BROWNIAN BRIDGE 11 − − − − − − − − − − − − − time in sec e rr o r γ = 0 . 45 (eq) γ = 0 . 25 (eq) γ = 0 . 45 (ad) γ = 0 . 25 (ad)order 1 / / Figure 3. Error given by (4) vs. average run-time for the strongapproximation of S ( w | B | ).5. Quantile Computation for a Weighted Reflecting BrownianBridge In this section we present an algorithm for the computation of the q -quantile of S ( w | B | ). The inputs of this algorithm, which will be denoted by Q ad ε ( w, q ), are theweight w , i.e., η and γ , as well as 0 < q < ε > 0. Thekey ingredient is the algorithm A ad n ( w, · ) from Section 4. Our construction is ratherad hoc and many improvements are possible. We mainly want to demonstratethe potential of using an adaptive algorithm as the building block for quantilecomputation.5.1. The Algorithm. The algorithm Q ad ε ( w, q ) starts with a precomputing stepin order to ideally determine the minimal integer n ∈ N such that e ( A ad n ( w, · )) ≤ ε. Due to the fast convergence of A ad n ( w, B ) towards the supremum S ( w | B | ), whichhas been observed in the numerical experiments in Section 4, we use E(∆ n ) with∆ n = | A ad2 n ( w, B ) − A ad n ( w, B ) | as an approximation to e ( A ad n ( w, · )). Moreover, we use a simple Monte Carlo al-gorithm X ( m ) n with m independent samples of ∆ n to approximate the expectationE(∆ n ). In the precomputing step we take m = 10 , and we compute the minimal integer n ∈ N of the form n = 10 · i with i ∈ N such that X ( m ) n ≤ ε. If no such n exists, then the output of the algorithm Q ad ε ( w, q ) is undefined. Weadd that this happens with probability zero if, as we conjecture, ∆ n p −→ k of independent samples of A ad n ( w, B ), which are independent from the precomputing step, too. The choice ofthis number is motivated by the following fact. Remark 12. Consider iid random variables Y , Y , . . . with a continuous densityfunction f . Assume that f ( F − ( q )) > q -quantile F − ( q ) of Y , and let Z q,k denote the d q · k e -th order statistic of Y , . . . , Y k . Then c q · √ k · ( Z q,k − F − ( q )) d −→ Z with Z ∼ N (0 , 1) and c q = f ( F − ( q )) p q · (1 − q ) , see, e.g., David and Nagaraja (2003, Thm. 10.3).Replacing the unknown constant c q from Remark 12 by one, we take k = d ε − e samples of A ad n ( w, B ), and the algorithm Q ad ε ( w, q ) returns the d q · k e -th orderstatistic of these samples.The following fact may be used to compute confidence intervals for the q -quantileof A ad n ( w, B ), which yields a quality control for the second step. Observe, however,that the precomputing step does not yield a rigorous link to the q -quantile of S ( w | B | ). Remark 13. Consider iid random variables Y , . . . , Y k , and let Y ( i ) the corre-sponding i -th order statistic. Furthermore, let 0 < α < 1, and assume that a, b ∈ { , . . . , k } with a < b satisfy P ( { Z ∈ { a, . . . , b − }} ) ≥ − α, where Z is binomially distributed with parameters k and q . Then [ Y ( a ) , Y ( b ) ] isa (conservative) confidence interval for the q -quantile of Y with confidence level1 − α , see, e.g., David and Nagaraja (2003, Sec. 7.1).5.2. Numerical Experiments. We compare the algorithm Q ad ε ( w, q ) and an al-gorithm Q eq ε ( w, q ) that is constructed as Q ad ε ( w, q ), but instead of A ad n ( w, · ) thenon-adaptive algorithm A eq n ( w, · ) is used as the building block. Furthermore, in theprecomputing step of Q eq ε ( w, q ) we fully use the findings from Figures 2 and 3 forfree in order to determine the value of n . This comes very close to choosing exactlyand at no computational cost the minimal integer n ∈ N such that e ( A eq n ( w, · )) ≤ ε ,and thus is very much in favor of Q eq ε ( w, q ) compared to of Q ad ε ( w, q ).We consider the error e ( Q ε ( w, q )) = E (cid:0)(cid:12)(cid:12) F − ( w, q ) − Q ε ( w, q ) (cid:12)(cid:12)(cid:1) (5)for Q ε ( w, q ) = Q ad ε ( w, q ) and Q ε ( w, q ) = Q eq ε ( w, q ), where F − ( w, q ) denotes the q -quantile of S ( w | B | ).We proceed as in the previous section. The only difference is that a deterministicquantity, F − ( w, q ), instead of a random variable, S ( w | B | ), is unknown in thedefinition of the error.The error e ( Q ε ( w, q )) and the average run-time are determined via a MonteCarlo simulation, where we use 10 replications for each value of ε , and the resultsare presented together with asymptotic confidence intervals as before. We also usethe same set of parameters η and γ , the same hardware system, and the sameprogramming language as before. The value of q is chosen as q = 0 . DAPTIVE QUANTILE COMPUTATION FOR BROWNIAN BRIDGE 13 In the Monte Carlo simulation we replace F − ( w, q ) by a highly accurate ap-proximation, namely 2 . γ = 0 . 25 and 2 . γ = 0 . 45. The latter isobtained as the d q · k e -th order statistic of k samples of A ad n ( w, B ), where n = 10 and where k is close to 4 · . According to the findings from Section 4 the dis-tributions of A ad n ( w, B ) and S ( w | B | ) should almost coincide for this value of n .Furthermore, using Remark 13 with confidence level 0 . 99, we have obtained theconfidence intervals [2 . , . γ = 0 . 25 and [2 . , . γ = 0 . γ .In Figure 4 we relate the actual error to the error tolerance ε . We see that bothalgorithms almost perfectly achieve the computational goal, namely, an error closeto the input ε . 10 − − − − − − ε e rr o r γ = 0 . 45 (eq) γ = 0 . 25 (eq) γ = 0 . 45 (ad) γ = 0 . 25 (ad) Figure 4. Error given by (5) vs. error tolerance ε for the quantileapproximation of S ( w | B | ) with q = 0 . Q ad ε ( w, q ) we observea polynomial order of convergence of about 1 / 2, while the corresponding order for Q eq ε ( w, q ) is only about 1 / 4. This corresponds to the findings from Section 4 andto Remark 12: The orders of convergence of A ad n ( w, · ) and A eq n ( w, · ) for the strongapproximation of S ( w | B | ) are given by ∞ and 1 / 2, respectively, see Figure 3,and the order of convergence for the quantile approximation should be 1 / 2, seeRemark 12.The algorithm Q ad ε ( w, q ) achieves an error 10 − in 5 seconds for γ = 0 . 25 andin 12 seconds for γ = 0 . 45, and even an error 10 − in less than 15 or 35 minutes,respectively. The algorithm Q eq ε ( w, q ) achieves the error 10 − in 2 minutes for γ =0 . 25 and (based on a extrapolation beyond the range of our numerical experiment)in more than 2 hours for γ = 0 . 45; the corresponding run-times for the error 10 − are 6 days and 6 years, respectively. − − − − − time in sec e rr o r γ = 0 . 45 (eq) γ = 0 . 25 (eq) γ = 0 . 45 (ad) γ = 0 . 25 (ad)order 1 / / / / Figure 5. Error given by (5) vs. average run-time for the quantileapproximation of S ( w | B | ) with q = 0 . Acknowledgement. The authors are grateful to Fabian Schlechter for program-ming and computer support. The authors thank Jan Henrik Fitschen, YannickKreis, and Lukas Mayer for valuable discussions. Stefanie Schwaar was supportedby the Deutsche Forschungsgemeinschaft (DFG) within the RTG 1932 ‘StochasticModels for Innovations in the Engineering Sciences’. References F. Akashi, H. Dette, and Y. Liu. Change-point detection in autoregressive modelswith no moment assumptions. Journal of Time Series Analysis , 39(5):763–786,2018.S. Asmussen, P. Glynn, and J. Pitman. Discretization error in simulation of one-dimensional reflecting Brownian motion. Ann. Appl. Probab. , 5(4):875–896, 1995.A. Aue and L. Horváth. Structural breaks in time series. J. Time Series Anal. , 34(1):1–16, 2013.A. Aue, G. Rice, and O. Sönmez. Detecting and dating structural breaks in func-tional data without dimension reduction. J. R. Statist. Soc. B , 80(3):509–529,2018.P. Bickel and M. Rosenblatt. On some global measures of the deviations of densityfunction estimates. Ann. Statist. , 1(3):1071–1091, 1973.A. Borodin and P. Salminen. Handbook of Brownian Motion . Springer, 2 nd edition,2002.J. M. Calvin. Average performance of a class of adaptive algorithms for globaloptimization. Ann. Appl. Probab. , 7:711–730, 1997.J. M. Calvin. A one-dimensional optimization algorithm and its convergence rateunder the Wiener measure. J. Complexity , 17(2):306–344, 2001.J. M. Calvin. Lower bound on complexity of optimization of continuous functions. J. Complexity , 20(5):773–795, 2004. DAPTIVE QUANTILE COMPUTATION FOR BROWNIAN BRIDGE 15 J. M. Calvin, M. Hefter, and A. Herzwurm. Adaptive approximation of the mini-mum of Brownian motion. J. Complexity , 39:17–37, 2017.T. H. Cormen, C. E. Leiserson, R. L. Rivest, and C. Stein. Introduction to Algo-rithms . MIT Press, 4 th edition, 2003.M. Csörgő and L. Horváth. Nonparametric methods for changepoint problems. In Quality Control and Reliability , volume 7 of Handbook of Statistics , pages 403–425. Elsevier, 1988.M. Csörgő and L. Horváth. Limit Theorems in Change-Point Analysis . Wiley seriesin probability and statistics. John Wiley & Sons Ltd., 1997.H. A. David and H. N. Nagaraja. Order Statistics . Wiley Series in Probability andStatistics. Wiley-Interscience, Hoboken, NJ, third edition, 2003.D. M. DeLong. Crossing probabilities for a square root boundary by a Besselprocess. Commun. Statist.-Theor. Meth. , 10(21):2197–2213, 1981.B. J. Eastwood and V. R. Eastwood. Tabulating weighted sup-norm functionalsof Brownian bridges via Monte Carlo simulation. In Asymptotic Methods inProbability and Statistics , pages 707–719. North-Holland, 1998.R. Fried and M. Imhoff. On the online detection of monotonic trends in time series. Biometrical Journal , 46(1):90–102, 2004.C. Gallagher, R. Lund, and M. Robbins. Changepoint detection in climate timeseries with long-term trends. Journal of Climate , 26(14):4994–5006, 2013.W. Härdle. Applied Nonparametric Regression . Cambridge University Press, Cam-bridge, 1990.M. Hefter and A. Herzwurm. Optimal strong approximation of the one-dimensionalsquared Bessel process. Commun. Math. Sci. , 15(8):2121–2141, 2017.Y. Hochberg. A sharper Bonferroni procedure for multiple tests of significance. Biometrika , 80(4):800–802, 1988.L. Horváth and P. Kokoszka. Inference for Functional Data with Applications .Springer, Berlin-Heidelberg-New York, 2012.C. Kirch. Resampling Methods for the Change Analysis of Dependent Data . PhDthesis, Universität zu Köln, Juni 2006.C. Kirch and J. Tadjuidje Kamgaing. Detection of change points in discrete-valuedtime series. In Handbook of Discrete-valued Time Series , pages 219–244. CRCPress, 2016.M. Orasch and W. Pouliot. Tabulating weighted sup-norm functionals used inchange-point analysis. J. Stat. Comput. Simul. , 74(4):249–276, 2004.E. S. Page. On problems in which a change in a parameter occurs at an unknownpoint. Biometrika , 44(1/2):248–252, 1957.K. Ritter. Approximation and optimization on the Wiener space. J. Complexity , 6(4):337–364, 1990.P. Salminen and M. Yor. On hitting times of affine boundaries by reflecting Brow-nian motion and Bessel processes. Period. Math. Hungar. , 62(1):75–101, 2011.S. Schwaar. Asymptotics for change-point tests and change-point estimators . PhDthesis, Technische Universität Kaiserslautern, 2016.S. Schwaar. A data-driven change-point estimator. arXiv:2010.12449 , 2020.M. Yor. On square-root boundaries for Bessel processes, and pole-seeking Brownianmotion. In Stochastic Analysis and Applications , volume 1095 of Lecture Notesin Math. , pages 100–107. Springer, 1984. Fachbereich Mathematik, TU Kaiserslautern, Postfach 3049, 67653 Kaiserslautern,Germany Current address : Fraunhofer ITWM, Fraunhofer-Platz 1, 67663 Kaiserslautern, Germany Email address : [email protected] Fachbereich Mathematik, TU Kaiserslautern, Postfach 3049, 67653 Kaiserslautern,Germany Email address : [email protected] R+V Versicherung AG, Raiffeisenplatz 2, 65189 Wiesbaden, Germany Email address : [email protected] Fachbereich Mathematik, TU Kaiserslautern, Postfach 3049, 67653 Kaiserslautern,Germany Email address : [email protected] Fraunhofer ITWM, Fraunhofer-Platz 1, 67663 Kaiserslautern, Germany Email address ::