On the connection between orthant probabilities and the first passage time problem
aa r X i v : . [ s t a t . C O ] J a n On the connection between orthant probabilitiesand the first passage time problem
E. Di Nardo
Dipartimento di Matematica, Universit`a degli Studi della BasilicataE-mail: [email protected]
Abstract
This article describes a new Monte Carlo method for the evalua-tion of the orthant probabilities by sampling first passage times of anon-singular Gaussian discrete time-series across an absorbing bound-ary. This procedure makes use of a simulation of several time-seriessample paths, aiming to record their first crossing instants. Thus, thecomputation of the orthant probabilities is traced back to the accuratesimulation of a non-singular Gaussian discrete-time series. Moreover, ifthe simulation is also efficient, this method is shown to be more speedythan the others proposed in the literature. As example, we make useof the Davies-Harte algorithm in the evaluation of the orthant prob-abilities associated to the ARFIMA(0 , d,
0) model. Test results arepresented that compare this method with currently available software.
Keywords:
First passage time, orthant probabilities, simulation of Gaus-sian discrete-time series, Davies-Harte algorithm.
For many computational problems in statistics and especially in economics(cf. for example Ku and Seneta, 1994), we have to compute the so-calledorthant probabilities (Tong, 1990) of the form P (cid:16) ∩ ki =1 { X i < S i } (cid:17) = Z O k π ) k/ | Σ | / exp (cid:26) − x T Σ − x (cid:27) d x , (1)where the random vector X ≡ ( X , X , . . . , X k ) follows a multivariate nor-mal distribution with zero mean and unity variance, Σ is a k × k symmetricpositive definite correlation matrix and O k = { x ≡ ( x , x , . . . , x k ) ∈ R k : x i < S i , S i ∈ R , i = 1 , , ..., k } k = 10 . Different numerical algorithms can be simply obtained via standard numeri-cal integrations, such as the classical Monte-Carlo method or some adaptivequadrature formulae. For k = 2 and k = 3 Genz (2001, see also referencestherein) sets up a competitive algorithm based on Gaussian quadratureswith adaptive integration. If k > , the infinite integration limits in (1)need to be carefully handled, either by using some type of transformationinto a finite region (see Genz, 1992 and 1993 and references therein) or byusing selected cutoff values. Genz has shown that the quadrature formulaetake less time and are more accurate than the Monte Carlo algorithm forsmall k, whereas the Monte Carlo method appears more efficient and hascomparable accuracy with k greater than 6 or 7 . In the last few years, a com-pletely different approach has spread, based on probability simulation. Ina review paper, Hajivassiliou et al. (1996) analyze the properties of severalsimulators and find that the one they propose with the label GHK (Geweke-Hajivassiliou-Keane) works better than all other methods by keeping a goodbalance between accuracy and computational costs. Here the idea consistsof estimating (1) by means of recursive conditioned probabilities involvingsamples of truncated standard normal random variables (r.v.’s). When Σhas a tridiagonal form or sparse structures, there are different numericalmethods available based on special decompositions of Σ (for similar matterssee Tong, 1990). When the orthant regions have bounds like S i = 0 for any i, the interest has been focused on expressing (1) in terms of the correlationcoefficients. Some solutions are available in closed forms for k = 2 and k = 3while for k ≥ k -dimensional recursion formula is found in Zhongrenand Kedem (1999) by using a polar coordinates trasformation.In short, with up to k = 10, there are different numerical methods avail-able by which orthant probabilities can be robustly and reliably computed2t low to moderate accuracy levels. High accuracy or high dimension prob-lems can require long computation times and it is still not clear what is thebest method for this type of problem.The aim of this paper is to explore the connection between the orthantprobabilities (1) and the first passage time (FPT) problem for a non-singularGaussian discrete time-series. This connection suggests a new way to evalu-ate the orthant probabilities by numerical methods. The keystone is choos-ing among some simulation procedures devoted to the construction of time-series paths. The mathematical kernel of the FPT problem for a non-singularGaussian time-series is recalled in section 2 (for a more detailed analysis seeDi Nardo, 2002). The Genz algorithm and the GHK simulator are summa-rized in section 3 together with some classical simulation procedures of thetime-series. As a working example, in the last section the orthant probabili-ties are computed in connection with the ARFIMA(0 , d,
0) model by meansof the FPT approach. Comparisons are also given with the results reachedvia the Genz algorithm and the GHK simulator.
Let us define the FPT r.v. of a discrete time-series X t through an absorbingtime dependent boundary S t as T = min t ∈ N { t : X t ≥ S t } . (2)In the following, suppose P ( X = x ) = 1 with x < S and in order tosimplify set x = 0 . First we observe that P ( T = k ) = P ( X < S , X < S , . . . , X k − < S k − , X k ≥ S k ) . (3)When X t is a Gaussian time-series with zero mean, unity variance and sym-metric positive definite correlation matrix Σ , it is P ( T = k ) = Z D k π ) k/ | Σ | / exp (cid:26) − x T Σ − x (cid:27) d x where D k = { x ≡ ( x , x , . . . , x k ) ∈ R k : x i < S i , i = 1 , , ..., k − , x k ≥ S k } . S ≡ ( S , S , . . . , S k ) and P k ( S , Σ) = P (cid:16) ∩ ki =1 { X i < S i } (cid:17) , equation(3) could be rewritten as P ( T = k ) = ( − P ( S , Σ) , if k = 1 ,P k − ( S , Σ) − P k ( S , Σ) , if k > P k ( S , Σ) = 1 − P ( T ≤ k ) . (4)The above equation provides the link between the orthant probabilities andthe FPT distribution of a non singular Gaussian time-series. This suggeststhe use of the machinery of the FPT problems for a speedy computing of(1).The next result allows us to estimate the FPT distribution function viasimulation of the time-series X t . Theorem 2.1
For an upper-bounded boundary S t , the FPT r.v. (2) is fair,i.e. P ( T < ∞ ) = 1 . Proof.
Denote by ρ ij the correlation coefficients of X i and X j for i, j =1 , , . . . , k. Set S max = max t ≥ S t < ∞ and ρ max = max i = j ρ ij . By using theSlepian inequality (Tong, 1990), i) if ρ max > P k ( S , Σ) ≤ P k ( S , Σ ρ max ) ≤ P k ( S max , Σ ρ max ) (5)where Σ ρ max is the matrix with 1 on the principal diagonal and ρ max elsewhere and S max is a vector whose k components are equal to S max . The random vector with normal distribution function N ( , Σ ρ max ) hasexchangeable r.v.’s and so (Tong, 1990) P k ( S max , Σ ρ max ) = Z ∞−∞ Φ k (cid:18) S max + √ ρ max z √ − ρ max (cid:19) φ ( z ) dz, S ∈ R where φ ( z ) = 1 √ π exp ( − z ) , Φ( x ) = Z x −∞ φ ( z ) dz, x ∈ R ; (6) ii) if ρ max ≤ P k ( S , Σ) ≤ P k ( S , I ) ≤ P k ( S max , I ) = Φ k ( S max ) (7)where I is the identity matrix and Φ( x ) is given in (6).4rom (5) and (7), it is lim k →∞ P k ( S , Σ) = 0 (8)and the result follows taking the limit in (4) for k going to infinity. Let us recall that the analytical results on FPT problems are mostly centeredon stochastic processes of diffusion type, where the Markov property playsa leading role in handling the related transition probability density function(pdf), see Ricciardi et al. (1999). The FPT distribution also has an explicitanalytical expression for the class of stochastic processes named the Levytype anomalous diffusion, in which the mean square displacement of thediffusive variable X t scales with time as t γ with 0 < γ < et al. , 2000).A typical simulation procedure lies in sampling N S values of the FPTr.v. T, by a suitable construction of as many time-discrete paths, and thento record the instants when such realizations first cross the boundary.Note that k consecutive observations of a standardized Gaussian time-series can be generated via the innovations algorithm from k i.i.d. standardnormal r.v.’s with O ( k ) operations (Brockwell and Davies, 1991). How-ever, if we add few hypotheses on the correlation function, faster methodsare available. For example the Durbin-Levinson algorithm (Brockwell andDavies, 1991) is an efficient procedure with O ( k ) operations if X t is sta-tionary and ρ t vanishing when t goes to infinity. If in addition the discreteFourier transform of the circulant embedding vector { ρ , ρ , . . . , ρ k − , ρ k , ρ k − , . . . , ρ } (9)is a positive real sequence, a still faster method is the Davies-Harte algorithm(Davies et al. , 1987) with computational cost O ( k log k ) (for a review onthe methods generating realizations of a Gaussian stationary process seePercival, 1992). 5herefore in order to evaluate the orthant probabilities (1), we suggestthe following Monte Carlo method: choose a faster method simulating thepaths of X t , record their first crossing instants trough the boundary S andthen estimate (1) trough (4).In the next section, a FORTRAN 77 implementation of this procedureis compared with the Genz algorithm and the GHK simulator that appearto be the more widespread methods computing orthant probabilities.The Genz method (Genz, 1993) transforms the original integral (1) intoan integral over unit hypercube P k ( S , Σ) = e Z e Z . . . Z e k Z d w , (10)where d w ≡ ( dw , . . . , dw k ) ,e i = Φ (cid:26) S c (cid:27) , for i = 1Φ ( S i − P i − j =1 c ij Φ ( − ( w j e j ) c ii ) , for i = 2 , · · · , k,c ij are elements of the lower triangular Cholesky decomposition of Σ , Φ( x )is given in (6) and Φ ( − ( x ) is its inverse. The idea of the Genz method is toapply a standard Monte Carlo method to (10) and evaluate { e i } by samplingpseudo-random numbers over (0 , . Let us observe that the Cholesky decom-position of Σ needs O ( k ) computations for a fixed k. The evaluation of theintegrand function in (10) takes O ( k k = 100 . The GHK simulator (Hajivassiliou et al. , 1996) estimates (1) by meansof ˜ P k ( S , Σ) = 1 N N X n =1 P ( A ) P ( A | y ,n ) · · · P ( A k | y ,n ; y ,n ; . . . ; y k − ,n ) (11)where { y i,n } are drawn sequentially from independent standard normal dis-tributions truncated to the events A i A i = ( −∞ < Y i < S i − P i − j =1 c ji Y j c ii ) c ij elements of the lower triangular Cholesky decomposition of Σ and { Y i } i ∈ N a sequence of i.i.d.r.v.’s having standard normal distribution. Againthe Cholesky decomposition of Σ takes O ( k ) computations for a fixed k, whereas the evaluation of (11) is of order O ( N × k ) for N sampling. Themethod has been implemented in FORTRAN 77 and in GAUSS, and theroutines are available at http://econ.lse.ac.uk/staff/vassilis/ up to k = 40 . The Davies-Harte algorithm generates a path with k +1 steps of a stationaryGaussian time-series having zero mean and autocovariances ρ , . . . , ρ k , suchthat the finite Fourier transform of (9) is non-negative. Craigmile (2003)has shown that the Davies-Harte algorithm is the most efficient way tosimulate an ARFIMA(0 , d,
0) process (Hosking, 1981). Let us recall thatthe standardized ARFIMA(0 , d,
0) time-series is Gaussian, stationary withzero mean and correlation function ρ k = d ( d + 1) · · · ( d + k − − d )(2 − d ) · · · ( k − d ) , k = 1 , , . . . (12)if ρ = 1 . To simulate a realization X , X , . . . , X k with the Davies-Hartealgorithm, the steps are: i) for n = 0 , , , . . . , k − g n = k − X j =0 γ ( j ) exp (cid:18) πijn k − (cid:19) + k − X j = k − γ (2 k − − j ) exp (cid:18) πijn k − (cid:19) ; (13) ii) for j = 0 , , . . . , k − X j +1 = 12 √ k − k − X n =0 Z n √ g n exp (cid:18) πijn k − (cid:19) (14)where Z , Z k − are real normal r.v.’s with zero mean and variance2 , { Z n } k − n =1 is a sequence of independent complex normal r.v.’s withindependent real and imaginary parts, each of variance 1 , and Z n = Z k − − n for n = k, . . . , k − . Indeed, inverting equation (13) it is ρ j = 12 k − k − X n =0 g n exp (cid:18) πijn k − (cid:19) , j = 0 , , . . . , k −
17o that
Cov ( X p +1 X q +1 ) = E ( X p +1 X q +1 )= 14( k − k − X n =0 2 k − X l =0 E ( Z n Z l ) √ g n g l exp (cid:18) πi [ pn − lq ]2 k − (cid:19) = 14( k − k − X n =0 q g k exp (cid:18) πi ( p − q ) k k − (cid:19) = ρ p − q . Then from (4), an estimation of (1) is˜ P k ( S , Σ) = 1 − number of first crossings before or equal kN S (15)where N S is the number of simulated paths.Observe that the step i) is evaluated just once for any fixed k withcomputational cost O ( k log k ) , whereas the step ii) is repeated N S timeseach with O ( k log k ) evaluations. To speed up the simulation with the fastFourier transform algorithm, round k up to the nearest power of two andtruncate the simulated series at the end.In the following, we take down some results in the evaluation of theorthant probabilities with the Genz algorithm, the GHK simulator and themethod suggested here. Specifically, Table 1 refers to the case of constantboundary S ( t ) = 1 , d = 0 . k ≥
20 and Table 2 refers to the case oflinear boundary S ( t ) = 2 − . t, d = 0 . k ≥ . In both cases, ˜ P k stands for numerical evaluations of (1).With regard to the Genz method, the maximum number N max of in-tegrand function evaluations has been allowed to k × , as suggested bythe same author. With this choice, the computed value ˜ P k has reached theabsolute accuracy 10 − required in input, except in some cases marked withstar in the tables. To the 99% confidence level, the estimated absolute errorhas been near 10 − for all k ≥ , except the cases marked with star in thetables where it has been of order 10 − . In Tables 1 ÷ CT G represents anevaluation in seconds of the time employed by the method in computing ˜ P k with such parameters.With regard to the GHK method, the number N of allowed simulationshas been set equal to N S as indicated in Tables 1 ÷ . In order to generatethe required N × k pseudo-random uniform numbers over (0 , , it has beenemployed the routine G05CAF of the NAG software library. In Tables 1 ÷ T G represents an evaluation in seconds of the time employed by the methodin computing ˜ P k with such parameters.The method here proposed has been implemented in a FORTRAN 77program on a server ALPHA with operating system OSF1 5.0 and by usingthe NAG software library: i) the routine G05FDF generates a vector of normal pseudo-random num-bers by means of the Box and Muller method (see for example Rubin-stein, 1981); ii) the routine C06HBF computes the discrete Fourier cosine transform ofthe sequence (13); iii) the routine C06EBF evaluates the discrete Fourier transform of the Her-mitian sequence in (14) with the fast Fourier transform (see Brigham,1974).The time window size has been set equal to 2 for k = 20 ÷
30 and equal to 2 for k = 31 ÷ . In Tables 1 ÷ , N S is the number of simulated sample pathswhereas CT F P T represents an evaluation in seconds of the time employedby the routine in computing ˜ P k with such parameters.As it would be theoretically expected, the GHK simulator and the Genzmethod cost more in computational time than the simulation procedure withthe Davies-Harte algorithm. Especially in the case of non-constant bound-ary, the Genz method is time-consuming if one would reach the requiredaccuracy. Furthermore, the proposed Monte Carlo method is particularlysuited to be implemented on a parallel computer in order to further cutthe computational time and to improve the approximation; this because thesample paths of the simulated time-series could be easily generated indepen-dently of one another. References [1] Balakrishnan V. (1985) Anomalous diffusion in one dimension.
PhysicaA , 569–580.[2] Brigham E.O. (1974)
The Fast Fourier Transform . Prentice-Hall.[3] Brockwell P.J., Davis R.A. (1991)
Time Series: Theory and Methods.
Springer-Verlag, New York, 2nd. ed.94] Craigmile P.F. (2003) Simulating a Class of Stationary Gaussian Pro-cesses Using the Davies-Harte Algorithm, with Application to LongMemory Processes.
Jour. Time Series Anal. , 505–511.[5] Davies R.B., Harte D.S. (1987) Test for Hurst effect. Biometrika ,95–101.[6] Di Nardo E. (2002) On first-passage problem for a non-singular Gaus-sian discrete-time series. Quad. Stat. , 51–70.[7] Di Nardo E., Nobile A.G., Pirozzi E., Ricciardi L.M., Rinaldi S. (2000)Simulation of Gaussian processes and first passage time densities eval-uation. Lecture Notes in Computer Science , 319-333.[8] Genz A. (1992) Numerical Computation of Multivariate Normal Prob-abilities.
J. Comp. Graph Stat. , 141–149.[9] Genz A. (1993) Comparison of Methods for the Computation of Mul-tivariate Normal Probabilities. Computing Science and Statistics ,400–405.[10] Genz A. (2001) Numerical Computation of Rectangular Bivariate andTrivariate Normal and t Probabilities. Preprint.[11] Gupta S.S. (1963) Probability Integrals of Multivariate Normal andMultivariate t. Ann. Math. Statist. , 792–828.[12] Hajivassiliou V., McFadden D., Ruud P. (1996) Simulation of Multivari-ate Normal Rectangle Probabilities and their Derivatives. Theoreticaland Computational results. J. Econom. (1-2), 85–134.[13] Hosking J.R.M. (1981) Fractional Differencing. Biometrika , 165–176.[14] Ku S., Seneta E. (1994) The number of peaks in a stationary sampleand orthant probabilities. Jour. Time Series Anal. , 385–403.[15] Molchan G.M. (1999) Maximum of a fractional Brownian motion: Prob-abilities of small values. Commun. Math. Phys. , 97–111.[16] Percival D.B. (1992) Simulating Gaussian Random Processes with Spec-ified Spectra.
Computing Science and Statistics , 534–538.[17] Plackett R.L. (1954) A Reduction Formula for Normal MultivariateIntegrals. Biometrika , 351–360.1018] Rangarajan G., Ding M. (2000) Anomalous diffusion and the first pas-sage time problem. Physical Review E , 120–133.[19] Rangarajan G., Ding M. (2000) First passage time distribution foranomalous diffusion. Physics Letters A , 322–330.[20] Ricciardi L.M., Di Crescenzo A., Giorno V., Nobile A.G. (1999) AnOutline of Theoretical and Algorithmic Approaches to First PassageTime Problems with Applications to Biological Modeling.
Math. Japon-ica , 247–322.[21] Rubinstein R.Y. (1981) Simulation and the Monte Carlo Method . Wiley,New York.[22] Sun H.J. (1988) A General Reduction Method for n -Variate NormalOrthant Probability. Comm. Stat. Theory and Methods , 3913–3921.[23] Tong Y.L. (1990) The Multivariate Normal Distribution . Springer-Verlag, New-York.[24] Zhongren N., Kedem B. (1999) On Normal Orthant Probabilities.
Chi-nese J. Appl. Probab. Statist. , 262–275.11able 1: Evaluations of the orthant probabilities ˜ P k for S ( t ) = 1 and d = 0 . . CT G , CT GHK and CT F PT represent an evaluation in seconds ofthe time employed in the computation of ˜ P k respectively for the Genzmethod, the GHK method and the FPT method. N S is the number ofthe simulated sample paths in the FPT method. Genz method GHK method FPT method k ˜ P k CT G ˜ P k CT GHK ˜ P k CT F P T N S
20 0.0924 2.91 0.0927 0.49 0.0925 0.38 200021 0.0835 7.19 0.0835 0.54 0.0838 0.41 210022 0.0756 5.11 0.0757 0.58 0.0752 0.41 210023 0.0684 8.07 0.0687 0.64 0.0686 0.42 220024 0.0620 3.69 0.0622 0.68 0.0618 0.42 220025 0.0563 8.98 0.0572 0.80 0.0558 0.46 240026 0.0511 9.45 0.0512 0.93 0.0513 0.50 265027 0.0463 4.28 0.0465 1.03 0.0460 0.50 265028 0.0422 6.81 0.0422 1.05 0.0423 0.50 265029 0.0383 7.18 0.0387 1.15 0.0385 0.52 275030 0.0349 4.93 0.0348 1.25 0.0335 0.53 280031 0.0317 3.33 0.0320 2.31 0.0318 1.43 390032 0.0289 5.37 0.0291 2.45 0.0288 1.47 395033 0.0264 5.57 0.0266 2.58 0.0262 1.48 400034 0.0240 8.95 0.0243 2.76 0.0232 1.48 400035 0.0220 6.03 0.0222 2.88 0.0215 1.48 400036 0.0199 1.07 0.0202 3.01 0.0198 1.48 400037 0.0182 6.64 0.0185 3.17 0.0182 2.30 630038 0.0167 4.37 0.0167 5.39 0.0165 2.30 630039 0.0153 6.95 0.0154 6.80 0.0155 2.34 640040 0.0140 3.09 0.0140 6.18 0.0137 2.38 650012able 2: