Effective persistency evaluation via exact excursion distributions for random processes and fields

Finding the probability that a stochastic system stays in a certain region of its state space over a specified time—a long-standing problem both in computational physics and in applied and theoretical mathematics—is approached through the extended and multivariate Rice formula. In principle, it applies to any smooth process multivariate both in argument and in value given that efficient numerical implementations of the high-dimensional integration are available. The computational method offers an exact integral representation yielding remarkably accurate results and provides an alternative method of computing persistency probability and exponent for a physical system. It can be viewed as an implementation of path integration for a smooth Gaussian process with an arbitrary covariance. Its high accuracy is due to efficient computation of expectations with respect to high-dimensional nearly singular Gaussian distributions. For Gaussian processes, the computations are effective and more precise than those based on the Rice series expansions and the independent interval approximation. For the benchmark diffusion process, it produces the persistency exponent that is essentially the same as the recently obtained analytical value and surpasses accuracy, interpretability as well as control of the error, previous methods including the independent or Markovian approximation. The method solves the two-step excursion dependence for a stationary differentiable Gaussian process, in both theoretical and numerical sense. The solution is based on exact expressions for the probability density for one and two successive excursion lengths. The numerical routine RIND computes the densities using recent advances in scientific computing and is easily accessible for a general covariance function, via a simple numerical interface. The work offers also analytical results that explain the effectiveness of the implemented methodology and elaborates its utilization for non-Gaussian processes.


Persistency of physical stochastic systems
The persistence phenomenon is fundamental in many areas of physical sciences. It is concerned with the fraction of time (or space) in which a physical system persists in a certain state. The persistence probability P(t) represents the probability that the system stays in the state in which it started at time zero. A remarkable fact about this probability is that in a wide range of physical contexts, it decays either algebraically P(t) ∼ t − θ , or exponentially P(t) ∼ e − θ t , the two being equivalent through the logarithm-exponent transformation, with θ in both the cases referred to as the persistence exponent. In a system in which the state space is continuous, persistency is typically referring to the system staying above a certain level. One then observes recurrences of the systems crossing this level and dividing the time (or space) into regions of excursions above and below the level.
It is not possible to fully account for the vast range of fields where the persistency plays an important role but to give a perspective on the importance of persistency determination we selected some of the recent contributions in various disciplines. In statistical physics, it has been used for non-equilibrium systems to characterize phase transitions, properties of the diffusion equation [1,2], survival of spin states [3], fluctuating steps, and interfaces [4]. In meteorology and climate research, persistency was used to statistically analyze the clustering of extreme weather events in historical climate records [5]. The recurrence interval distribution in earthquakes has brought persistency analysis into seismology [6] and analysis of heartbeat intervals demonstrates its importance in physiology [7]. In an analysis of stocks during economic crises and distress, the tail of the distribution interval between return volatilities above a certain threshold proves that the persistency coefficient is an important characteristic that allows assessing financial risks of investments [8,9]. Other financial questions appear in option pricing under Lévy random walk where crossing time distributions for upper or lower barriers are of interest [10]. All these examples highlight the importance of the distribution of recurrence times or return intervals between extreme events for characterization and understanding of the behavior of physical systems and phenomena.
For stochastic processes, the persistence probability can be viewed as the ratio of two path integrals. The numerator integrates over paths obeying the condition of persistence and the normalizing denominator integrates over all paths. The approach through path integrals has its limitation and in [4], we read 'But as yet there is no general scheme for an arbitrary correlation function for calculating the persistence exponent.' While the statement remains still valid if one considers analytical evaluations of the persistency, it has been largely overlooked that there is a computational method utilizing the generalized Rice formula that allows effective computations not only of the persistency exponent but also persistency probabilities. The method works for any covariance in a Gaussian process and extends to non-Gaussian models obtained by subjecting the Gaussian measure to a random time change. We advocate the method as a computational tool for studying persistency problems in physical systems.

The problem and some of its early history
The problem of finding the intensity of the zero crossings of a random function has been probably first formulated as a mathematical question by Littlewood and Offord [11]. About the same time it was solved by Rice [12], who presented the celebrated Rice formula for the first time, also pointing out the physical context: 'Although this problem is of some physical interest I have been unable to find references to any earlier work. Problems of this nature occur in the investigation of the current reflected by small random irregularities along telephone transmission lines.' The result is also known under the name Kac-Rice formula due to its solution for random polynomials in [13]. One of the central problems in Rice's second article on random noise [14] is the statistical characterization of the zeros of a stationary Gaussian process. Rice's formula for the expected number of zeros, and more generally, of non-zero level crossings, is the first step, but Rice also presents the in-and exclusion series, the 'Rice series', for the distribution of the time between two successive mean level crossings.
The distribution of the number of zero crossings is naturally connected to the distribution of the successive lengths of excursions above and below zero. Both problems were followed during the decades following Rice's article. However, elegant analytical solutions are elusive and the following quote by Rice [14] remains valid: 'The problem of determining the distribution function for the distance between two successive zero seems to be quite difficult and apparently nobody has as yet given a satisfactory solution.' The original Rice series was improved considerably by Longuet-Higgins who obtained a rapidly converging moment series for the probability density of zero-crossing intervals [15,16] and also compared approximations based on the initial terms in the series with experimental results [17] and with earlier alternative series, suggested by McFadden [18,19].
Some studies by McFadden and Rainal [19][20][21] are of particular interest for the present article since they present systematic theoretical as well as experimental studies of the dependence between successive crossing intervals. Three approximations were studied, independence, 'quasi'-independence, which assumes that the sum of two successive intervals is independent of the next one, and the Markov type dependence. The first two cases were analyzed by renewal type arguments using the Laplace transform [22,23], and some numerical solutions were compared to experiments. The Markov assumption [24] was tested by variance and correlation parameters against experiments. All three assumptions were rejected for general Gaussian processes.
On the theoretical side, formulas for crossing moments of arbitrary order under minimal assumptions for Gaussian processes have been obtained [25], and a generalized Rice formula for not necessarily Gaussian processes has been presented by Zähle [26]. Further, a regression technique was introduced in [27,28] as the first step towards the numerical algorithms that are discussed in this paper.

Renewed interest and new exact tools
In the last two decades, the interest in the tail behavior of excursion times remerged in material science, optics, statistical physics, and other areas [29]. The emphasis was on the tail distribution rate referred to as the persistency exponent. The 'independent interval assumption' (IIA) was applied both to Gaussian processes and to other process models and compared to experiments [23]. The diffusion processes were analyzed in general [1,30,31] and an analytical solution for the persistency coefficient for the diffusion of order two was obtained by a rather deep combination of arguments across different developments in theoretical physics [32].
During the same years, new tools were developed in applied probability and scientific computing, and the exact formula for the joint density of successive crossing intervals was given [33]. It involves the conditional expectations of the derivatives at crossings and the indicator functions that the process stays above or below the level between crossings. No analytic expressions are known, but the integral forms are readily computable by the high-dimensional integration made possible by a carefully designed numerical approach.

Persistency of non-equlibrium systems-the diffusion case
Nonequilibrium dynamics of physical systems is one area of physics where the persistency problem was investigated through empirical studies, for example, for first-passage exponent in two-dimensional soap froth [34], and in a liquid crystal system exhibiting Ising-like behavior [35], the propagation of slow-combustion fronts in paper sheets [36], for the temporal fluctuations of complex crystalline structures [37], for variation in the spin orientation of a sample of laser-polarized 129 Xe gas [38]. However, due to the statistical nature of the phenomena, experimental results are rare and the theoretical understanding of the persistency phenomena is critical. It is confirmed by theory and experiment that the relation between the persistency exponent and the characteristic defining the decay of temporal or spatial correlation functions is complex. Probably the most studied case is a diffusing field starting from a random initial configuration [31,38]. This model is attractive both theoretically due to the linear form of the diffusion equation in any dimension, table 1, and experimentally because of its natural physical interpretation in various non-equilibrium systems. For these reasons, the model constitutes a mathematically convenient and physically important benchmark for any methodology aiming at the persistency assessment. This benchmark is used throughout the paper and the results for it obtained through the proposed computational methodology are compared with the results existing in the literature , table 2. Examples of what can be expected for non-diffusive systems are also presented in the paper. We also present the potential of implementation of pathwise integrals for the persistency in non-Gaussian case that is important when non-linear generalizations of the diffusion equation are considered [36].
Type: Noise and sea waves

Theoretical foundation of the approach
To get a proper insight into the methodology, let us recall a generalized Rice formula. We have opted for a fairly general while mathematical formulation that is a minor generalization of Theorems 6.2 and 6.4 in [39]. The reader will see in Example 1 below and section 3 that this is actually a very practical result. Let Z: U → R d be a twice continuously differentiable Gaussian field, U ⊂ R k , k d, with a non-singular covariance matrix at any point. Then for the k − d dimensional measure  of the subset  u ( ) L of these tʼs in Λ ⊂ U such that Z(t) = u: In addition, assume that for each t ä U one has another continuous field, Y t : W → R n , so that (t, w) → (Z(t), Y t (w)) defined on U × W is Gaussian. Then, for a bounded and continuous function g: U × C(W, R n ) → R: When the mathematical formalism is stripped out, there are several key points to be made about equation (1) and, more generally, equation (2). The inner integration is essentially about averaging over the multivariate Gaussian distribution. Consequently, efficient evaluation of such averages is the central problem in making the formula usable. One does not need to restrict to the Gaussian universe and the results, under some technical restrictions, generalize to more general stochastic processes. However, the integration in the conditional expectation significantly adds to the computational challenge. In the Gaussian case, the integration challenges have been addressed by proper use of computational linear algebra and numerical integration techniques in [40,41], where a practically useful routine for computation of high-dimensional normal integrals has been developed. Brodtkorb [42] combined all the described ideas into a powerful tool, adding new tests to control the accuracy, and embedding it into a user-friendly code for Gaussian process crossing problems. The routine, called RIND, based on these methods was included in the MATLAB toolbox WAFO [43,44].
A WAFO-tutorial [45] guides the user through the technicalities of using the computational tools. It also gives information about the main computational engine, the routine RIND, and presents examples of the code. An illustration how the result can be practically used is presented next to study the backward and forward delays via Rice's formula.

Example 1. We derive the joint pdf of A B
, for the clipped process in figure 1. We consider the zero-level crossing case, but the general case can be treated similarly. In this problem,] Hence, by stationarity of X, the joint probability density of B A , is given by The persistence problem for non-Gaussian continuous-time models is extremely challenging and not much work has been done in this direction, see [46]. While, in general, path integrals obtained from the generalized Rice formula are still conceptually valid within the smooth but non-Gaussian domain, their efficient evaluation needs to be addressed. The applicability of the approach to non-Gaussian models has not been intensively studied but one specific case of the conditional Gaussian models has been addressed [47,48], namely Laplacetype non-Gaussian models built by Brownian motion with a random (Lévy) time change. A natural method here is to apply conditioning on the underlying random time and then do additional integration over it. Such an approach is computationally more expensive and has additional challenges due to the effect of conditioning on the crossing value. Nevertheless, it is still computationally manageable and in some special cases random time change analytical 'shortcuts' allows for significant improvement of the computational effectiveness. In the final section of this work, we provide a short overview of this technique.
2. Persistence of a stochastic process 2.1. The point process of crossing instants and the IIA The excursion time distribution can be formulated through the point process of crossing instants. We consider a smooth process X(t), − ∞ < t < ∞ , and the time instants of u-level crossings, S i , leading to two sequences of interlaced intervals of lengths T i + , i = ± 1, ± 2, ..., for the excursions above u-level and T i -, i = ± 1, ± 2, ..., for the analogous excursions below u. We label the interval that contains the origin [ − B, A], with a forward delay time A to the first crossing on the positive side, and a backward delay time B since the last crossing on the negative side, while its length T Figure 1 explains the principle for indexing. The variable δ is introduced to keep track of excursions above, δ = 1, or below u, δ = − 1. The process D c (t) = + 1( − 1) when X (t) > u( < u) is called the clipped version of the X-process at level u. If the smooth process X(t) is stationary, the clipped process D c (t) is also stationary. The point process of u-level crossings {S i }, is a stationary point process.
The clipped process is not easy to study due to the dependence between A, B, and T i ʼs. In an approach that is quite popular in computational physics [31] one approximates this point process by a simpler renewal process with independent crossing intervals. The covariance of the clipped Gaussian process is readily available and the covariance of the simple renewal process connects to the distribution of independent intervals. Matching the covariances of the two processes may give some insight into the distribution of the excursions times. Indeed, some accurate approximations of the exponent of the tail distribution of T i ʼs have been obtained by this method.
The match goes via the Laplace transforms and is based on the following relation between the covariance R of the renewal process with its Laplace transform R and the probability distribution Φ of independent intervals with its Laplace transform (in the probabilistic sense, i.e. the Laplace transform of a probabilistic measure) denoted by Ψ: where μ is the average length of the renewal interval. Connection is then made by plug-in-style substitution of the covariance function of a clipped process for the covariance R in the above. The above relation is in agreement with formula (215) in [31]. However, the success of this method seems to end just right there, as the approach does not provide any insight what are the probability distributions that satisfy these equations or even when they exist at all. It has some applications to study persistency exponent since exponential asymptotics can be observed through the largest negative pole of the Laplace transform. Any other questions about the excursion distribution would require not only inverting the Laplace transform but first evaluating the Laplace transform of the covariance. We conclude that the method has quite limited applications.

Dependence of level crossing intervals
The IIA is a popular method of approximating the persistence of stochastic processes, thus it is important to investigate the degree of the dependence between crossing intervals and evaluate its effect on the approximation accuracy. It is well-known that no Gaussian process can have exactly independent mean level crossing intervals [15,19,49]. In this sense, the IIA is never exact. Here, we give an account of numerical evaluation of the degree of the dependence for Gaussian processes. The exact approach through the generalized Rice formula is explained and evaluated using the RIND routine, the technical details of the computations are given in [45]. We compare the exact RIND pdf with the IIA-based pdf for the spectra in table 1, extending the cases extensively studied in [50], and illustrate the results. We group these examples in two groups: spectra with near independent halfperiods, i.e. intervals between two subsequent crossings, and with strongly dependent half-periods. Note that the marginal distributions are computed from the exact pdf by integration and not by the IIA approach, which is incapable of retrieving the marginal distribution as discussed later. The name convention for the spectra and covariances in table 1 is as follows. LHk and WHk hint at [15] and [51]. The diffusion spectra BMSd were used in [31]. WN is low-frequency white noise with a Butterworth approximation BS, and the Jonswap spectrum J is an example of an ocean wave spectrum. Note that the spectra and covariance functions are listed in the unnormalized form.
In the examples, they are normalized so that the variance of the process and its derivative are both equal to zero and thus the average length of a zero-crossing interval equals π.
We start with a class of processes where successive zero-crossing intervals are 'almost' independent, representing diffusion in d dimensions, where d = 1, 2 correspond to physically feasible experiments [31], section 9. For all quantile levels, except for the most extreme, the true pdf and the one obtained by the independent approximation are almost identical as shown in figure 2 Left. Similar results have been seen for other d-values.
The remaining graphs in figure 2 show joint pdf for spectra with near independent intervals. There is a clear visual agreement between the exact pdf (red) and the pdf with independent margins. The Gaussian spectrum WH0 and the approximating LH1 are centered at zero frequency and follow the IIA model almost perfectly. The white noise spectrum WN deviates noticeably from how one normally envisages independent variables. For the approximating Butterworth spectrum BS the pdf is near to the independence. Figure 3 Left-Middle shows examples with clear or even strong dependence. As seen in the graphs, the dependence can take many different shapes, which makes it difficult to catch it in a simple parametric form.
The RIND function is not limited to mean level crossings but works for any level as illustrated on the shifted Gaussian spectrum model WH6 for levels u = 0.1, 0.3. Figure 3 Right shows dependence and good agreement between the RIND results and simulated pdf:s.

The persistence via the exact distributions
The most common meaning of persistence is as the tail of the first-crossing distribution  processes with non-negative correlation function, while for oscillating correlation only upper and lower bounds have been found. Non-negative correlation:- [52], Thm. 1.6 gives a precise meaning to the 'exponential tail' property: if the correlation function r(t) of a stationary Gaussian process X(t) is everywhere non-negative, there exists a non- The limit b r is necessarily finite [53] and thus it is meaningful to formulate the persistence tail as The IIA was used in [1] to find θ IIA for different dimensions of the diffusion processes defined by the covariances in table 1, and compared with simulations, while experimental evidence for d = 1 was presented in [38]. An efficient simulation procedure to estimate the persistence for arbitrary dimension was devised in [55], which yielded the simulation-based evaluation of the corresponding persistence exponents. Recently, a definite answer for dimension d = 2 was given in [32], namely θ(2) = 3/16 = 0.1875, a value that agrees with the simulation in [55]. At present, no exact values are known for other dimensions.
The persistence estimates in [55], Tab. 1 come from a large and well-controlled simulation experiment of over 10 8 realizations of first crossing events. We have designed the following numeric procedure to estimate the persistence exponent.
(i) Fix a maximum T m and a grid t k = kΔ t , t m Δ t = T m ; this resembles the logarithmic grid in the [55] simulations. T m should be chosen so large that any numeric or stochastic uncertainty in the tail is filtered out.
(ii) Compute, with RIND, Q T for T = t k , k = 1, K, m; in [55] Q T is estimated by simulation.
(iii) Make repeated independent runs to compute Q T and take the average Q ; T RIND uses Monte Carlo integration for extreme cases, and averaging will reduce the stochastic uncertainty.
(iv) Make a regression of Q log 2 T ( ) against T; this can be done globally, over the entire interval [0, T m ], or locally to allow for slow asymptotics in equation (8). The persistence exponent is then estimated as the slope of the fit. We used quadratic polynomial fit to identify any typical trend in the exponential e ( θ+ A( T)) T .
We used this scheme to estimate the persistence exponent for the 10 dimensions in [55]. We set Δ t = 0.05 and T m = 15 for d 4, which corresponds to the time span used in that paper, Δ t = 0.1, and T m = 30, 30, 20 for d = 1, 2, 3. We computed Q T as the average of 400 independent runs of RIND with the SOBNIED method and highest precision. The quadratic fit indicated that the local slope increased slowly with T. Table 2 shows the exponents from [55], labeled NL and the RIND-values, represented by the local θ T in the middle of the interval, T = T m /2.
For dimensions d = 1, 2 we give two values in parenthesis to take account of the asymptotic character of the persistence exponent. To accomplish this, we increased T m to 40, 35, respectively, and estimated the local rate of decay for large T. For d = 1, the value θ = 0.1206 is the best estimate over the interval (0, 30), while the value 0.1203 is the stable value for large T. The NL-value 0.1205 is probably too high, based on a too short time span. The RIND-value 0.1203 agrees with the value computed by IIA. For d = 2 the stable value is 0.1875, equal to the theoretical value 3/16, given in [32]. For the diffusion in two dimensions, we evaluate the persistence for other than zero levels and obtain values 0.2052, 0.2283, 0.2426, for levels u = 0.1, 0.2, 0.3, respectively, suggesting some monotone dependence of the persistency on the level.

The Durbin-Rychlik formula
We have seen the IIA has limited application to studies of the excursion intervals even in the cases when the independence seems nearly satisfied. On the other hand, the exact method not only yields more accurate results but is also readily available. The root of the exact formula is Slepian's 'doubly conditioned' process [28], which explicitly describes a process with a level upcrossing at 0 and downcrossing at t; the conditioning is implicit in Rice's original paper [14]. The Rice series approximations achieve the qualifications by restricting the number of extra crossings in the interior of the intervals by higher-order moments for the number of crossings. Recursive formulas how to compute all truncated moments to a very high computational cost are found in [56]. The exact formula for the distribution of level crossing intervals in Gaussian processes, developed in [57,58], can be fairly easily derived from (2). In our derivations, we use the following notational conventions from [33]: X s,t,u = (X(s), X(t), X(u)), X X s X t X u s t u , , Moreover, a X s,t b means that for each u ä (s, t): a X(u) b.
To understand how the formula is used for the crossing distribution one has to reformulate the problem. We consider pairs of instants of an u-upcrossing and the following u-downcrossing and ask first how many of the upcrossings, on average, are found in the unit interval. This is the intensity of the upcrossings and thus is found from the standard Rice formula X X u f u r r u r where r X is the (smooth) covariance function of X. Finding the survival function of the excursion time reduces thus to the evaluation of how many of these pairs on average represent the length larger than, say, t 0 > 0. This is a more challenging task, which can be approached by embedding it into a higher dimension and applying the generalized Rice formula.
Namely, in the theorem, we set n = 2 and take Z(t) = Z(t, s) = (X t , X t+s ), Λ = [0, 1] × [0, t 0 ], Y t (w) = X(w), and g(t, Y t ) = {X t,t+s > u}. It should be pointed out that in this case the function g is not continuous and some technical results should be used to augment (2); see the proof of theorem 7.1 in [59] for mathematical details of the argument. We note that the measure is of the dimension zero and thus it is again the counting measure. Now, when it comes to Z t ( )  it is a 2 × 2 triangular matrix with X X , t t s ( )   + on the diagonal. Then using our abbreviated notation, we have Z X t det t t s , ( )   = + and it is not difficult to see that (2) yields the searched average where the second equality is by stationarity. Consequently, the pdf of T + is given by An analogous argument is given in the appendix to obtain the distribution of A, B. For a Gaussian process, when u is equal to the mean level, the clipped process is symmetric with respect to the abscissa. Let μ + and μ − denote the mean value of T i + ʼs and T i -ʼs, respectively, while μ is their common value in the symmetric case. For a symmetrically clipped process, the joint density of the forward delay A and the backward delay B has the density For the asymmetric case, with δ = ± 1 indicating the status of the interval that contains the origin, the distribution of (A, B, δ) is given through where f A,B|δ stands for the conditional density. The discrepancy between the long-run interval distributions in a stationary point process and distributions taken from a frozen starting point has been first discussed for telephone calls. The solution has been worked out using the Palm measures, renewal processes, horizontal window conditioning, and the Rice formula. The mathematical foundations have been long resolved [60][61][62][63], as in the key renewal theorem. We observe that the intervals of a constant value (plus or minus) that include the origin of the horizontal axis are not distributed the same as the similar intervals (excursions) as observed over the entire real line, see figure 4. As seen in the graphs, the excursion time distribution over the whole line is closer to the distribution of the distance from the origin to the first crossing rather than the distribution of the entire excursion interval containing the origin. Statistically speaking the origin of the horizontal line hits larger intervals than those following the excursion time distribution. Thus the well-known inspection paradox, according to which observing a renewal interval at a fixed time yields on average longer length of it than that of an average renewal interval, can be also demonstrated in the dependent renewal times case.
We note that the exact distribution can be dealt with in the same manner over other-than-zero levels at which the process is clipped. Using the method, we have explored the relationship between the persistency and the crossing level for the diffusion process in dimension two having the covariance r t t 2 sech 2 X ( ) ( ) = . In figure 5 (left), we see increasing dependence of the persistency on the crossing level.

The joint crossing distribution and Markovian approximation
The investigation of the excursions can be extended to the joint distribution of subsequent intervals. In particular, it can shed some light on the dependence structure and the accuracy of the IIA as well as of the Markovian approximation.
For our purpose, we add an extra condition to Slepian's 'doubly conditioned' process to get a 'triply conditioned' process with a zero downcrossing at t = 0 with upcrossings at s < 0 < u. We derive the exact joint distribution from (2) by the embedding method. In the notation of that theorem, t = (r, t, s), r < 0 < s, Like most, so-called, 'explicit solutions' to mathematical problems, the expectation in (19) has to be evaluated numerically. The complexity is the same as computing the distribution of the maximum of a smooth Gaussian process X(t), [ ]  P , see [64] for efficient software, and [65] for an analysis of numerical accuracy. Due to the strong local dependence for smooth Gaussian processes, the way to compute the expectation in To obtain sufficient accuracy one may have to take a dense grid which can result in an almost singular multivariate Gaussian distribution. The result, the MATLAB routine RIND, is included in the package WAFO [44]. The approach gives upper and lower bounds for the approximation error with controlled accuracy of computations.
We have seen that the relation between the dependence and the accuracy of the IIA approximation is complicated and rather sensitive to small departures from the independence assumption. In what follows, we have compared the exact RIND pdf with the IIA-based pdf for all the models in table 1 and we present numerical dependence measures. First, we report that for the diffusion spectra BMSd used in [31], we have nearly independence as can be seen in figure 2. In figure 5 (right), we present the dependence for shifted Gaussian spectra and in table 3 for other examples of spectra. The name convention for the spectra and covariances is as follows. LHk and WHk hint at [15] and [51]. WN is low-frequency white noise with a Butterworth approximation BS, and the Jonswap spectrum J is an example of an ocean wave spectrum.
For each of the spectra, we computed the correlation between successive half periods and the Kullback-Leibler distance between the exact pdf and the IIA pdf. Figure 5  We conclude this part by providing an insight into the Markovian approximation (MA) for zero crossing intervals, which dates back to McFadden [19] and Rainal [20,21], who also tested the model using the sequence of correlations. Having the exact joint distribution of successive zero crossing intervals a natural next step is to test the Markov chain dependence of the full sequence of crossing intervals. If f t t , T T , 1 2 1 2 ( ) is the joint density of two successive crossing intervals, one can construct a Markov transition kernel as k t t f t t f t , Since the numerical algorithm gives the joint density in the discretized form it is natural to construct a discrete Markov chain with discrete states x 1 , x 2 , K, x n, and transition matrix T x We consider a Markov process with transition probabilities given by equation (20) and use it to obtain an approximation for the whole crossing sequence. A natural, exact, and rather strong test of the model can be based on the trivariate distribution of three consecutive intervals. The exact tri-variate density f t t t , , 1 2 3 ( ) is computed by RIND with the same degree of complexity as the bivariate density. If the crossing sequence is Markovian, then the conditional pdf of the length T 2 , given the length T 1 of the previous interval, should be equal to its conditional pdf, given the lengths of the two previous ones (T 0 , T 1 ), i.e., for all We illustrate the technique on the WH3, WH6, WH8 processes, as examples of processes with 'almost independent', 'moderately dependent', and 'strongly dependent' successive intervals. We use the RIND applications cov2ttpdf and cov2tttpdf to compute the 2D and 3D densities in equations (20) and (21) and then compute and plot the left-hand sides conditional densities for selected values of x j and, for each x j , for different x i . Figure 6 shows the results. From the graphs, we can conclude the dependence between intervals for WH3. However, the colored curves in each subplot do not deviate much from the respective black curve, except for the shortest interval, which indicates that the Markov chain approximation can be used. For WH6 the deviation from a Markov model is stronger, and for WH8 the Markov model fails altogether.

Evaluation of the excursions in higher dimensions
The generalized Rice formula is not limited to one dimension and can be applied to random fields depending on temporal and multivariate spatial arguments. In higher dimensions, persistence is harder to define since the excursion sets cease to have a simple form of a union of compact intervals and may create complex random topological structures. Figure 7 presents excursion sets for fields W(t, p) with a time variable t and a space variable p. Although both fields may have similar behavior along the temporal axis, the excursion sets and persistence may be quite different if the spatial relations are taken into account. For example, if the intensity of the crossings is measured by the average length of the contour per unit square, it is quite clear that the second field has longer contours than the first one. Despite the apparent complexity, some properties can be still quantified by the generalized Rice formula and RIND routine.
As an illustration, let us consider the point process on the planes at the local maxima of the process. This problem was thoroughly studied in [66]. In further development, the stochastic representation of the process at the local maximum known as the Slepian model has been given in [67]. Here we consider the local maxima within the excursion sets together with the closest point at the crossing contour at which the tangent line is perpendicular to the contour. The disk located at the maximum and tangent to the contour at that point is entirely enclosed in the excursion set. If such disks were used for the renormalized field to obtain the same crossing rates in space and time, it would lead to ellipsoids and their sizes quantify the persistence of the stochastic field, see figure 8 (left). In [68], the method of the generalized Rice formula and the RIND integration has been applied to evaluate the three-dimensional distribution (A, R, Θ), where A is the height in the local maximum, R is the distance to the point at the level contour, and Θ is the spatial-temporal azimuth of the line connecting the location of the maximum with the point on the contour. The distribution is measured by the observed relative frequencies of the local maxima with a given property.  A is the value of the field at the maximum, R is the arrow's length, and Θ is the angle the arrow makes with the temporal axis. We see the joint distributions of (R, A) (middle) and (Θ, R) (right) given that A is bigger than the 83%-percentile of its distribution.
We proceed in the familiar fashion by first addressing the overall frequency of the positive local maxima obtained from equation (1) Marking additionally by the positive value of the local maxima allows evaluating the intensity of the positive local maxima per temporal-spatial unit The algebraic formula in terms of special functions has been obtained through a rather lengthy derivation in [66], pp. 351-353. However, to obtain the joint density of (R, A, Θ) the problem is considerably more complex and there is no hope for an analytical solution. Instead, one can resort to the RIND methodology as follows.
We consider W(q), which plays the role of Y t in (2). The event that marks the positive local maxima is then , where  r is the disk centered at 0 having radius r. In agreement with our previous notational conventions, let B be also identified with an indicator function of an event B, i.e. equal one if B holds and zero otherwise. Then, by embedding into dimension four with where W θθ , W r are partial derivatives in the polar coordinates t p r r , c o s , s i n .
The derivation is similar to the cases already presented and thus omitted.
For the density f (h, r, θ) in equation (24), RIND is utilized and an example of computations for the JONSWAP spectrum spread by a cosine function are shown in figure 8 (middle-right). There we see the distribution of rather rare events due to the condition on A to be above the 83% quantile. The RIND computed distribution matches the intense simulation data well.
It follows from the asymptotic properties of the Gaussian fields that as the value of A is getting larger, the evaluation of the joint pdf of (A, R) ) is getting easier (and faster). This is because, in the numerical evaluation, the indicator B  in the formula can be replaced by one. This is in contrast to using the simulation approach to approximate the pdfs involving A by counting observed frequencies-very high local maxima occur very seldom and the computational cost increases.

Non-Gaussian stochastic models
Although the generalized Rice formula is also valid for non-Gaussian processes, the RIND is developed strictly for Gaussian models and the computational extensions beyond the Gaussian domain are challenging even for one-dimensional excursion sets. However, there is one class of models for which RIND is still applicable although the computational challenges are greater than in the Gaussian case. In these models, the Brownian noise has the time affected by a random time change represented by a non-decreasing Lévy process. Technically speaking, this process, given the random time change, becomes a non-homogeneous Gaussian process to which both (2) and RIND apply. Thus by introducing additional integration with respect to random time one can still utilize the convenience of high-dimensional strongly dependent multivariate Gaussian distributions. Things are not that straightforward though, since one has to address first what is the distribution of the random time change around a crossing of the process. This was considered in [47] and a stochastic form of the noise at the crossing instants was given. At the moment, the numerical implementation is not yet available in a convenient interface to address the excursion distribution problem. Instead, we illustrate the difference between Gaussian and non-Gaussian models featuring Lévy time change through a simpler case of the slope distribution at the crossing. The problem does not involve high-dimensional integration and is reduced to finding the joint distribution of the process and its derivative at a crossing.
It is well known that, for a stationary Gaussian process, the derivative at a level crossing is independent of the level value, which may seem counterintuitive. One of the few analytical results in the crossing distribution theory states that the derivative at a crossing is Rayleigh distributed. These two elegant results are no longer true for non-Gaussian models. To illustrate this we consider a non-Gaussian response to the Brownian motion with time changed by the gamma process. In figure 9, we show the derivative distributions at different crossing levels for processes having the same covariance of a smoothed Ornstein-Uhlenbeck type, which are obtained from the kernels (impulse responses) presented in the graphs as well. The computations for the non-Gaussian case have been made using the Rice formula and numerically evaluating the joint distribution of the process and its derivative at the single location [48,69]. The derivative becomes steeper at higher levels, which is due to the steepness of the chosen kernel around zero and thus kernels of different shapes will produce different results.

Conclusions
The progress in solving analytically the classical problem of the excursion distribution has been slow. In computational physics, some quasi-analytical approaches, such as the IIA or Markovian approximation of the persistence exponent, have been developed. However, their applicability is limited to the properties that are easily expressible through the Laplace transforms On the other hand, recent advances in statistical computing have made it possible to compute probabilities and expectations for very high-dimensional while nearly singular Gaussian distributions. It has been discussed how the MATLAB implementation RIND of these methods can be used to solve intricate level crossing problems for stochastic processes and fields. The approach is very precise for the excursion sets of Gaussian processes in one dimension. Moreover, as shown through examples from dynamical random fields and certain non-Gaussian models, a wide range of advanced excursion problems are now numerically accessible. Consequently, the Rice-based method combined with effective high-dimensional integration routines constitutes as an important computational tool for studying the level crossing distributions of a large variety of random physical systems.