{'fraction_non_alphanumeric': 0.0585131150085364, 'fraction_numerical': 0.012222567127114697, 'mean_word_length': 4.142258579409417, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Large-Deviation Functions for Nonlinear Functionals of a Gaussian Stationary Markov Process

We introduce a general method, based on a mapping onto quantum mechanics, for investigating the large-T limit of the distribution P(r,T) of the nonlinear functional r[V] = (1/T)\int_0^T dT' V[X(T')], where V(X) is an arbitrary function of the stationary Gaussian Markov process X(T). For T tending to infinity at fixed r we find that P(r,T) behaves as exp[-theta(r) T], where theta(r) is a large deviation function. We present explicit results for a number of special cases, including the case V(X) = X \theta(X) which is related to the cooling and the heating degree days relevant to weather derivatives.


I. INTRODUCTION
The "persistence" of a continuous stochastic process has been the subject of considerable recent interest among both theoreticians and experimentalists in the field of nonequilibrium processes. Persistence is the probability P (t) that a stochastic variable x(t) of zero mean does not change sign up to time t. Systems studied include reaction-diffusion processes, phase-ordering kinetics, fluctuating interfaces and simple diffusion from random initial conditions [1]. Experimental measurements have been carried out on breath figures [2], liquid crystals [3], soap froths [4], and diffusion of Xe gas in one dimension [5]. In "coarsening" systems like these, which do not possess a definite length or time scale, the persistence has a power-law decay, P (t) ∼ t −θ , at late times, and the persistence exponent θ is in general non-trivial. In such systems, the normalized two-time correlation function, C(t 1 , t 2 ) = x(t 1 )x(t 2 ) /[ x 2 (t 1 ) x 2 (t 2 ) ] 1/2 has the "scaling" form, C(t 1 , t 2 ) = f (t 1 /t 2 ), depending only on the ratio of the two times. In such systems, a simplification is achieved by introducing the logarithmic time-scale T = ln t, and the normalized variable X(T ) = x(t)/ x 2 (t) 1/2 , since the correlation function of X depends only on the time difference T 2 − T 1 , i.e. the process X(T ) is stationary. Thus one is led to consider stationary stochastic processes. These processes are, of course, also of interest in their own right. In the new time variable, the persistence decays exponentially, P ∼ exp(−θT ).
In this paper we consider the integrated quantity where V (X) is an arbitrary function of the stochastic variable X(T ). Hence r is a functional of V . For the special case V (X) = H(X), where H(X) is the Heaviside step function, r is just the fraction of time for which X(T ′ ) > 0 in the time interval O ≤ T ′ ≤ T . In this case, the probability distribution, P (r, T ), of r for given T is just the "residence -time" distribution which, together with the related "sign-time" distribution, where V (X) = sgn(X), has attracted a lot of attention earlier in the mathematics literature [6,7] and more recently in the statistical physics community [8]. Here we consider a general function V (X). We restrict our attention, however, to the class of processes where X(T ) is a stationary Gaussian Markov process, for which some analytic progress can be made. A case of some interest, which we analyze in detail, is V (X) = XH(X). This case is relevant to weather derivatives as we now explain. It has been estimated that about 1 trillion dollars of the 7 trillion dollar US economy is weather related [9]. For example, weather conditions directly affect agricultural outputs and the demand for energy products and indirectly affect retail businesses [10]. Weather derivatives, first introduced in 1997, are financial instruments which allow hedging (by, for example, energy suppliers) against adverse weather conditions over a period of time.
Here "adverse" might mean an unusually warm winter, when low demand for energy would affect the supplier's profits, as well as an unusually cold one when the supplier is unable to meet the demand. Temperature derivatives, the most common form of weather derivatives, are based on the concepts of "heating degree days" and "cooling degree days", which are (rough) measures of the cumulative demand for heating and cooling respectively.
Let X(T ) be the temperature at time T in a given city. On a given day, n, the mean, X n , of the highest and lowest temperatures is recorded. The number of cooling degree days (CDD), over a period of N days, is given by CDD = N n=1 max(X n − X o , 0), where X o is a reference, or baseline, temperature, while the number of heating degree days (HDD) is HDD = N n=1 max(X o − X n , 0). In the present paper we will, for simplicity, use an integral over continuous time rather than a sum, so that the cool-ing degree days are given by CDD = where H(X) is the Heaviside step function. Thus CDD is the integrated temperature excess (over the reference temperature) restricted to those periods where the temperature is above the reference level. It is a crude measure of the amount of cooling (air conditioning) required during the period T and also of the energy required to produce this cooling. Note that the power consumption of an air conditioner actually varies, for small temperature differences, as the square of the temperature difference between the room and ambient temperatures, so a better measure of the energy required would be In the notation of this paper, (1). The number of heating degree days (HDD) is similarly given by , which is a measure the amount of heating required in the period T and of the energy required to produce this heating. Estimating the likelihood of large deviations from the mean in CDD or HDD is clearly of interest to energy companies.
To make further progress, a realistic statistical model of temperature fluctuations is required. At this stage, however, any realistic model is too intractable to allow analytical progress. To illustrate the general approach we will instead employ a simple, though unrealistic, model in which the temperature, X(T ), is taken to be a stationary Gaussian Markov process. We will discuss the limitations of this model, and possible improvements, in the Conclusion. As a further simplification, we take the reference temperature X o to be the mean of X(T ), though this simplification is not essential and can be relaxed.
The outline of the paper is as follows. In section II we introduce the general approach to the problem of computing the distribution, P T (r), of the quantity r defined by Eq. (1). Using a path-integral representation, the calculation of P T (r) is mapped onto a problem in quantum mechanics, in which the function V (X) appears in the potential energy. In the limit of large T , only the quantum ground state energy is required. The final result takes the form P T (r) ∼ exp[−θ(r)T ], where the function θ(r) is a large-deviation function. For the case of the signtime distribution, corresponding to V (X) = sgn(X), r lies in the range −1 ≤ r ≤ 1, and it is clear that θ(1) is just the usual persistence exponent, θ, since r = 1 requires X(T ′ ) > 0 for 0 ≤ T ′ ≤ T . In section III, the method is illustrated on a number of special cases, of which V (X) = αX and V (X) = βX 2 /2 are studied first as exactly-soluble tutorial illustrations before turning to the CDD problem (in the simple form outlined above) and finally the sign-time distribution. The last two examples can be solved analytically in various regimes, and numerically elsewhere. Section IV concludes with a discussion and summary of the results.

II. GENERAL APPROACH
Consider the general stationary Gaussian Markov process (Ornstein-Uhlenbeck process) dx/dt = −µx + ξ(t), where ξ(t) is Gaussian white noise with zero mean, and correlator ξ(t)ξ(t ′ ) = 2D δ(t − t ′ ). After the change of variables t = T /µ, x = (2D/µ) 1/2 X, the equation takes the form where η(T ) = 0 and The probability distribution of X(T ′ ) for 0 ≤ T ′ ≤ T is given by whereẊ(T ′ ) ≡ dX/dT ′ , and N is a normalization constant. Our goal is to calculate the probability distribution of r, defined by equation (1). In practice it is convenient to look at the distribution P u (u) of the quantity u = rT . Its Laplace transform is where Z(s) is given by the path integral We are interested in the limit T → ∞. It is convenient to impose periodic boundary conditions, X(T ) = X(0), since this restriction will not change the results in the large-T limit. Furthermore, the exponential in (4) should strictly contain the combination (Ẋ + X) 2 instead of (Ẋ 2 +X 2 ). The missing term, 2XẊ, is however a perfect derivative, whose integral vanishes for periodic boundary conditions. Finally, with these boundary conditions the function Z(s) is the imaginary-time Feynman pathintegral that gives the partition function of a quantum particle with Hamiltonian H = p 2 /2 + X 2 /2 + sV (X) at inverse temperature T , p being the canonical momentum conjugate to X. For T → ∞ the ground state dominates to give, in this limit, where E(s) is the ground-state energy for the Schrödinger equation with potential For s = 0 the problem reduces to a simple harmonic oscillator, and E(0) = 1/2. The stochastic process x(t) studied above corresponds to the position of a Brownian particle in an external potential µx 2 /2. For the case of a pure Brownian motion (µ = 0), Kac derived a formalism [6] to study the distributions of arbitrary functionals V [x] which also used a mapping to the Schrödinger equation. Note, however, that the method presented above for the µ = 0 case differs in details from the original Kac formalism.
To illustrate the method we discuss two simple examples, before turning to some nontrivial cases, including the CDD problem.
For this case we have This case is actually trivial since u, being a sum of zeromean Gaussian variables, is itself a zero-mean Gaussian variable. All we require, therefore, is the variance, given by For the Ornstein-Uhlenbeck process (2), with noise correlator (3), one easily finds X(T 1 ) X(T 2 ) = (1/2) exp(−|T 1 − T 2 |). Inserting this result in (11), and extracting the leading large-T behavior, gives u 2 → α 2 T , and therefore r 2 → α 2 /T . Hence the asymptotic distribution of r is given by (neglecting prefactors) We now show how the general machinery we have set up in section II recovers this result. The potential U (X) in the Schrödinger equation (8) takes the form This is just a harmonic oscillator with a shifted origin, so E(s) = 1/2 − s 2 α 2 /2 and, using (7), for large T . To recover P T (r) we can invert the Laplace transform as follows. Neglecting pre-exponential factors, This integral can, of course, be evaluated exactly. Here, however, we use the method of steepest descents, which is valid for large T and can be readily generalized to the other, less trivial, cases that we will discuss. Writing the integrand in the form exp[T g(s)], we have g(s) = rs + s 2 α 2 /2. The integral is dominated by the saddle point at s = −r/α 2 , where g = −r 2 /2α 2 . The integration contour is deformed to pass over the saddle point, which lies on the real s axis. The saddle point is thus a minimum of g(s) with respect to variations of s along the real axis. The final result, ignoring non-exponential prefactors, is identical to (12). The structure of (12) is with We can easily generalize this method to arbitrary V (X). The path-integral approach gives the general result (neglecting prefactors) where Using the steepest-descent method for the integral gives where we have inserted E(0) = 1/2. The next example is a simple application of this idea.
This case illustrates the power of the method. The potential energy is now so we have a harmonic oscillator again but with a modified frequency, ω = (1 + βs) 1/2 , giving E(s) = (1 + βs) 1/2 /2. Thus θ(r) = max Note that θ(r) now has its minimum at r = β/4, which is just the mean value of r [noting that X 2 = 1/2 follows from (2) and (3)], while large (r → ∞) and small (r → 0) values of r are strongly suppressed. An expansion of θ near its minimum value gives a Gaussian distribution P T (r) ∼ exp[−4T (r − β/4) 2 /β 2 ]. This agrees with the result expected from the central limit theorem, i.e. P T (r) ∼ exp[−(r − µ) 2 /2σ 2 ], with mean µ = β/4 and variance σ 2 = β 2 /8T . For T → ∞, the distribution becomes very narrow such that, at fixed r, the central limit theorem fails to give the correct asymptotics. The form (16) thus gives the behavior in the extreme tails of the distribution at large T . The function θ(r) is a "largedeviation function" which controls the distribution of r for large T . We note that this special case with V (X) = βX 2 /2 was also studied recently by Farago [11] in the context of power fluctuations in the Langevin equation (2) by a somewhat different method. The probability density function of the dissipated power in reference [11] is precisely the distribution P T (r) studied here and the corresponding large-deviation function θ(r) has the same expression as in Eq. (22). However our derivation seems simpler and easily generalizable to other forms of V (X) as we show below.
C. V (X) = XH(X): "Cooling Degree Days" In this case the quantity rT gives the integrated value of X over the interval T , restricted to those values where X > 0. If X is the excess temperature over some baseline value where cooling becomes necessary, rT is a crude measure of the total energy consumption required to provide the cooling [V (X) = X 2 H(X) would be a better measure, as discussed in the Introduction]. In the present simple model we have taken the mean temperature equal to the reference temperature (with X(T ) being the deviation from the mean), though this is not an essential restriction. The Schrödinger equation (8) takes the form where ψ ′′ ≡ d 2 ψ/dX 2 . The required solutions can be expressed in terms of parabolic cylinder functions, D p (x), using the standard solutions of the parabolic cylinder equation y ′′ −(x 2 /4+a)y = 0 [12]. Selecting the solutions that satisfy the physical boundary condition ψ(±∞) = 0 gives where A, B are constants and The ratio A/B and the energy eigenvalues E(s) are obtained from matching the wavefunction and its derivative at X = 0, i.e. we require ψ + (0) = ψ − (0) and ψ + ′ (0) = ψ − ′ (0). This yields the eigenvalue equation The determination of E(s) from these equations is not possible analytically for general s. However, in the regimes s → ∞, and s → −∞, which determine the small-r and large-r behavior respectively of θ(r), analytical progress is possible. We consider these in turn.

The limit r → ∞
As we shall see, in this limit it is sufficient to compute E(s) in the limit s → −∞. This can either be done directly from (29), or using the following (simpler) physical argument.
Recall that, for V (X) = XH(X), the potential U (X) in the Schrödinger equation is For s → −∞, the potential has a deep minimum, of depth −s 2 /2, located at X = −s. The wavefunction is exponentially small at X = 0, ψ(0) ∼ exp(−s 2 /2), and the change in the potential in the regime X < 0 has an exponentially small effect on the ground state energy. Thus The same result can be obtained, after some algebra, directly from (29). Inserting the result in (20) gives θ(r) = max s (−s 2 /2 − rs) = r 2 /2. Since the maximum occurs at s = −r, the calculation is self-consistent for r → ∞. Thus we obtain This large-r result has the same form as equation (17), which gives the general-r result for potential V (X) = αX (with α = 1 here). This correspondence is not so surprising: For r → ∞, the dominant processes X(T ′ ), for 0 ≤ T ′ ≤ T , will correspond to large positive X, so the step function H(X) in V (X) plays no role in this limit.

The limit r → 0
We shall see that this limit corresponds to s → ∞. Again, one can use physical arguments as a short cut. For s → ∞, one has U (X) = X 2 /2 for X < 0, with essentially a hard wall at X = 0. This gives, to leading order, E(s) → 3/2. We need, however, the leading correction to this result. Since the wave function does not penetrate far into the wall, we can neglect the X 2 term in the potential for X > 0, i.e. write U (X) = sX. The Schrödinger equation then simplifies to The wavefunction for X ≤ 0 is again a parabolic cylinder function, while for X ≥ 0 it can be expressed as an Airy function: where A, B are constants, and p − = E(s) − 1/2 as before. Matching the wavefunction and its first derivative at X = 0 gives the eigenvalue equation where Γ(x) is the Gamma function. In the limit s → ∞, the left-hand side of (37) tends to infinity, so the righthand side (RHS) must also diverge in this limit. The ground state corresponds to the first divergence, where p − → 1. Therefore we write p − = 1 − ǫ in (37), and seek the leading behavior as ǫ → 0. This gives RHS → −(2/π) 1/2 /ǫ. Inserting this result in (37) gives, to leading order, Note that the value of s at which the maximum occurs in (39) is s = (a/3r) 3/4 , justifying our use of a large-s analysis of E(s) for the limit r → 0.
That θ(r) → 1 for r → 0 is intuitively clear, since r = 0 requires X(T ) ≤ 0 for all T . This reduces to the usual persistence probability of the Markov process (2), for which θ = 1.
3. θ(r) for r near r Equations (32) and (39) give analytical results for θ(r) in the limits of large and small r respectively. For general r, θ has to be computed numerically. There is, however, one other regime where analytical progress is possible, namely for r close to its mean value, where the central limit theorem (CLT) applies.
The mean value is given by The Gaussian distribution for X gives immediately X = 0 and |X| = 2/π X 2 1/2 = 1/ √ π. Thus In a similar way, the variance of r can be calculated, by exploiting the Gaussian properties of the process X(T ). A tedious but straightforward calculation gives, for T → ∞, The central limit theorem then gives the behavior of P T (r) for r near r c as P T (r) = ( Inserting the explicit expressions for r and correct to leading nontrivial order in (r − r ). The full result for θ(r) can be obtained by numerically solving (29) for the ground state energy E(s) for each value of r, then using (20) to find the corresponding θ(r). The result is displayed in Figure 1, with the asymptotic forms for r → ∞ and r → 0 indicated. Note the very sharp rise to the value unity as r → 0, as indicated by Eq. (39). In terms of the CDD problem, the behavior below the minimum (i.e. for r < r ) determines the probability of an unusually small CDD ("cool summer"), while the region above the minimum corresponds to an unusually large CDD ("hot summer"). The fact that the function θ(r) initially increases less rapidly to the right of the minimum than to the left indicates that (within this very simple model) summers with a slightly larger than average CDD are more probable than those with a slightly smaller than average CDD. This asymmetry is a consequence of the nonlinear relation between CDD and the temperature fluctuations, which are symmetric about the mean in our model. It should be stressed that the integration period T has been taken to be large, to justify the steepest descent calculation. In practice this means that T (the length of a summer, say) must be large compared to the correlation time of the temperature (a few days, perhaps), which is not totally unreasonable. The structure of the calculation is similar to that of the preceding subsection. The Schrödinger equation is The solutions are parabolic cylinder functions, where now Matching the wavefunction and its derivative at X = 0 gives the eigenvalue equation which, using standard identities relating D p (0) and D p ′ (0) to gamma functions [12], reduces to Although this equation cannot be solved analytically for general s, the limits s → 0 and s → ∞ are tractable.
Note that E(s) = E(−s) by symmetry, so there is no need to consider s → −∞ separately. The analysis starts from the potential well, U (X) = X 2 /2 + s sgn(X). For small s, the ground state energy E(s) is perturbatively close to, but slightly smaller than, E(0) = 1/2. Therefore we write where we anticipate δ(s) = O(s 2 ) from the symmetry E(s) = E(−s). Inserting this form in (49), (51) becomes Expanding to fourth order in s and second order in δ gives The maximum occurs at s = −r/2 ln 2 + O(r 3 ), so our study at small s is self-consistent at small r. For s → ∞, on the other hand, the potential develops a hard wall at the origin, and has a depth of −s next to the wall. Therefore we write E(s) = −s + 3/2 − ǫ, with ǫ small. Putting this form in (49), (51) becomes Taking the limits s ≫ 1, and ǫ ≪ 1 readily leads to to leading order for large s. Using this in (20) gives The maximum occurs at s = [2 √ π(1 + r)] −2/3 ], which tends to infinity as r → −1 so the calculation is selfconsistent in this limit. The symmetry of the problem under r → −r leads to the more general result The function θ(r) is plotted in Figure 2, with only the region r ≥ 0 shown explicitly. The limiting behavior for small r and r close to unity is also shown.

IV. CONCLUSION
In this paper we have presented a general method for computing the asymptotic behavior of the the probability distribution of the quantity r = (1/T ) where V (X) is an arbitrary function and X(T ) is an Ornstein-Uhlenbeck stochastic process representing the motion of a Brownian particle in the presence of a stable parabolic potential µX 2 /2. The main new result is that for µ > 0, the distribution of r, for large window size T has the form P T (r) ∼ exp[−θ(r)T ], where θ(r) is the large-deviation function. The calculation proceeds via a mapping onto a quantum mechanical problem described by the Schrödinger equation (8) for a particle moving in the potential (9), where s is the Laplace variable conjugate to rT . The inverse Laplace transform can be performed using steepest descents in the limit T → ∞.
Two non-trivial applications have been presented. The first is a calculation of the large-deviation function θ(r) for the cooling degree days problem, based on a simple model of temperature fluctuations. The model used assumes that the temperature fluctuation X(T ) is described by the Ornstein-Uhlenbeck process (2), i.e. that X(T ) is a stationary Gaussian Markov process with timeindependent noise strength. This may not be a realistic model for several reasons: (i) A simple Markov process is not thought to be an optimized model of temperature fluctuations, which tend to exhibit stronger autocorrelations than a Markov process. A more plausible statistical model, according to [10], writes the fluctuation, X n , of the mean (average of daily high and low) temperature on day n as the linear combination X n = k m=1 w n−m X m + η n , where the "memory kernel" w n−m is a decreasing function of n − m and η n is uncorrelated Gaussian noise. The Markov case corresponds to k = 1. Here X n is the difference between the measured mean temperature on day n and its expected value. The latter should contain a 365 day seasonal variation (roughly sinusoidal) (ii) The noise strength should also contain a 365 day seasonal variation: the variance of the temperature fluctuations can be different at different times of the year. (iii) The reference temperature for the calculation of CDD and HDD should in general be different from the expected temperature. We hope to incorporate some of these features in future studies.
The second non-trivial application is to the "sign-time" distribution. In the context of "sign-time", we point out that the asymptotic form of the "sign-time" distribution has a markedly different behavior in the Ornstein-Uhlenbeck process (where a Brownian particle moves in a stable parabolic potential µX 2 /2) studied here compared to the ordinary Brownian motion (µ = 0). In the later case, the "sign-time" distribution P T (r) is independent of the window size T for all T and is given by P T (r) = 1/π √ 1 − r 2 [7]. In the former case (µ > 0), on the other hand, the "sign-time" distribution depends explicitly on the window size T and P T (r) ∼ exp[−θ(r)T ] for large T where the large-deviation function θ(r) has been computed exactly in this paper. We further note that for this "sign-time" problem, the function θ(r) can also be obtained using the "independent interval approximation" (IIA) [13,14], which exploits the fact that the intervals between zero crossings are statistically independent for a Markov process. In fact, for renewal type processes where the successive intervals are statistically independent, the "sign-time" distribution has been computed by Godrèche and Luck [15] using the interval size distribution as an input. Their result can be simply extended to calculate the "sign-time' distribution for the Ornstein-Uhlenbeck process. The IIA also has the virtue that it can be used to obtain approximate results for non-Markov processes. Persistence exponents, for example, are often given rather accurately by the IIA [1]. However, it is not straightforward to adapt this IIA method to general nonlinear functions V (X), whereas the path-integral approach and mapping onto quantum mechanics adopted here is readily applicable to any V (X).
Here we have only considered Gaussian Markov processes mainly because they are simple and amenable to analytical calculations. Recently the calculations of the asymptotic distributions for the "sign-time" and other related quantities such as "local-time" have been extended to non-Gaussian Markov processes where a Brownian particle moves in an arbitrary stable or unstable potential and moreover exact results have been obtained [16] when the underlying potential is random as in the Sinai model. The extension of these methods and results presented here to non-Markov processes, however, still remains as one of the outstanding challenges for the future.
V. ACKNOWLEDGMENT SM thanks A. Comtet for useful discussions.