Networks of ·/G/∞\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\cdot /G/\infty $$\end{document} queues with shot-noise-driven arrival intensities

We study infinite-server queues in which the arrival process is a Cox process (or doubly stochastic Poisson process), of which the arrival rate is given by a shot-noise process. A shot-noise rate emerges naturally in cases where the arrival rate tends to exhibit sudden increases (or shots) at random epochs, after which the rate is inclined to revert to lower values. Exponential decay of the shot noise is assumed, so that the queueing systems are amenable to analysis. In particular, we perform transient analysis on the number of jobs in the queue jointly with the value of the driving shot-noise process. Additionally, we derive heavy-traffic asymptotics for the number of jobs in the system by using a linear scaling of the shot intensity. First we focus on a one-dimensional setting in which there is a single infinite-server queue, which we then extend to a network setting.


Introduction
In the queueing literature, one has traditionally studied queues with Poisson input. The Poisson assumption typically facilitates explicit analysis, but it does not always align well with actual data; see, for example, [11] and references therein. More specifically, statistical studies show that in many practical situations, Poisson processes underestimate the variability of the queue's input stream. This observation has motivated research on queues fed by arrival processes that better capture the burstiness observed in practice.
The extent to which burstiness takes place can be measured by the dispersion index, i.e., the ratio of the variance of the number of arrivals in a given interval, and the corresponding expected value. In arrival streams that display burstiness, the dispersion index is larger than unity (as opposed to Poisson processes, for which it is equal to unity), a phenomenon that is usually referred to as overdispersion. It is desirable that the arrival process of the queueing model takes the observed overdispersion into account. One way to achieve this is to make use of Cox processes, which are Poisson processes, conditional on the stochastic time-dependent intensity. It is an immediate consequence of the law of total variance that Cox processes do have a dispersion index larger than unity. Therefore, this class of processes makes for a good candidate to model overdispersed input processes.
In this paper we contribute to the development of queueing models fed by input streams that exhibit overdispersion. We analyze infinite-server queues driven by a particular Cox process, in which the rate is a (stochastic) shot-noise process.
The shot-noise process we use is one in which there are only upward jumps (or shots) that arrive according to a homogeneous Poisson process. Furthermore, we employ an exponential 'response' or 'decay' function, which encodes how quickly the process will decline after a jump. In this case, the shot-noise process is a Markov process; see [16, p. 393]. There are several variations on shot-noise processes; see, for example, [10] for a comprehensive overview.
It is not a novel idea to use a shot-noise process as stochastic intensity. For instance, in insurance mathematics, the authors of [5] use a shot-noise-driven Cox process to model the claim count. They assume that disasters happen according to a Poisson process, and each disaster can induce a cluster of arriving claims. The disaster corresponds to a shot upwards in the claim intensity. As time passes, the claim intensity process decreases, as more and more claims are settled.
Another example of shot-noise arrival processes is found in the famous paper [13], where it is used to model the occurrences of earthquakes. The arrival process considered in [13] has one crucial difference with the one used in this paper: it makes use of Hawkes processes [9], which do have a shot-noise structure, but have the special feature that they are self-exciting. More specifically, in Hawkes processes, an arrival induces a shot in the arrival rate, whereas in our shot-noise-driven Cox model these shots are merely exogenous. The Hawkes process is less tractable than the shot-noisedriven Cox process. A very recent effort to analyze ·/G/∞ queues that are driven by a Hawkes process has been made in [8], where a functional central limit theorem is derived for the number of jobs in the system. In this model, obtaining explicit results (in a non-asymptotic setting), as we are able to do in the shot-noise-driven Cox variant, is still an open problem.
In order to successfully implement a theoretical model, it is crucial to have methods to estimate its parameters from data. The shot-noise-driven Cox process is attractive since it has this property. Statistical methods that filter the unobservable intensity process, based on Markov Chain Monte Carlo (MCMC) techniques, have been developed; see [3] and references therein. By filtering, they refer to the estimation of the intensity process in a given time interval, given a realized arrival process. Subsequently, given this 'filtered path' of the intensity process, the parameters of the shot-noise process can be estimated by a Monte Carlo version of the expectation maximization (EM) method. Furthermore, the shot-noise-driven Cox process can also be easily simulated; see, for example, the thinning procedure described in [12].
In this paper we study networks of infinite-server queues with shot-noise-driven Cox input. We assume that the service times at a given node are i.i.d. samples from a general distribution. The output of a queue is routed to a next queue, or leaves the network. Infinite-server queues have the inherent advantage that jobs do not interfere with one another, which considerably simplifies the analysis. Furthermore, infinite-server systems are frequently used to produce approximations for corresponding finite-server systems. In the network setting, we can model queueing systems that are driven by correlated shot-noise arrival processes. With regards to applications, such a system could, for example, represent the call centers of a fire department and police department in the same town.
The contributions and organization of this paper are as follows. In this paper we derive exact and asymptotic results. The main result of the exact analysis is Theorem 4.6, where we find the joint Laplace transform of the numbers of jobs in the queues of a feedforward network, jointly with the shot-noise-driven arrival rates. We build up toward this result as follows. In Sect. 2 we introduce notation and we state the important Lemma 2.1 that we repeatedly rely on. Then we derive exact results for the single infinite-server queue with a shot-noise arrival rate in Sect. 3.1. Subsequently, in Sect. 3.2, we show that after an appropriate scaling the number of jobs in the system satisfies a functional central limit theorem (Theorem 3.4); the limiting process is an Ornstein-Uhlenbeck (OU) process driven by a superposition of a Brownian motion and an integrated OU process. We then extend the theory to a network setting in Sect. 4. Before we consider full-blown networks, we first consider a tandem system consisting of an arbitrary number of infinite-server queues in Sect. 4.1. Then it is argued in Sect. 4.2 that a feedforward network can be seen as a number of tandem queues in parallel. We analyze two different ways in which dependency can enter the system through the arrival process. Firstly, in Model (M1), parallel service facilities are driven by a multidimensional shot-noise process in which the shots are simultaneous (which includes the possibility that all shot-noise processes are equal). Secondly, in Model (M2), we assume that there is one shot-noise arrival intensity that generates simultaneous arrivals in all tandems. In Sect. 5 we finish with some concluding remarks.

Notation and preliminaries
Let ( , F, {F t } t≥0 , P) be a probability space, in which the filtration {F t } t≥0 is such that (·) is adapted to it. A shot-noise process is a process that has random jumps at Poisson epochs, and a deterministic 'response' or 'decay' function, which governs the behavior of the process. See [16,Section 8.7] for a brief account of shot-noise processes. The shot noise that we use in this paper has the following representation: where the B i ≥ 0 are i.i.d. shots from a general distribution, the decay function is exponential with rate r > 0, P B is a homogeneous Poisson process with rate ν, and the epochs of the shots, that arrived before time t, are labeled t 1 , t 2 , . . . , t P B (t) . As explained in the introduction, the shot-noise process serves as a stochastic arrival rate to a queueing system. It is straightforward to simulate a shot-noise process; for an illustration of a sample path, consider Fig. 1. Using the thinning method for nonhomogeneous Poisson processes [12], and using the sample path of Fig. 1 as the arrival rate, one can generate a corresponding sample path for the arrival process, as is displayed in Fig. 2. Typically, most arrivals occur shortly after peaks in the shot-noise process in Fig. 1, as expected.
We write (i.e., without argument) for a random variable with distribution equal to that of lim t→∞ (t). We now present well-known transient and stationary moments of the shot-noise process, see Appendix 2 and, for example [16]: with B distributed as B 1 , We remark that, for convenience, we throughout assume (0) = 0. The results can be readily extended to the case in which (0) is a nonnegative random variable, at the cost of a somewhat more cumbersome notation.
In the one-dimensional case, we denote β(s) = E e −s B , and in the multidimensional case, where s = (s 1 , s 2 , . . . , s d ), for some integer d ≥ 2, now denotes a vector, we write The following lemma will be important for the derivation of the joint transform of (t) and the number of jobs in system, in both single-and multi-node cases.
Proof See Appendix 1.

A single infinite-server queue
In this section we study the M S /G/∞ queue. This is a single infinite-server queue, of which the arrival process is a Cox process driven by the shot-noise process (·), as defined in Sect. 2. First we derive exact results in Sect. 3.1, where we find the joint transform of the number of jobs in the system and the shot-noise rate, and derive expressions for the expected value and variance. Subsequently, in Sect. 3.2, we derive a functional central limit theorem for this model.

Exact analysis
We let J i be the service requirement of the i-th job, where J 1 , J 2 , . . . are assumed to be i.i.d.; in the sequel J denotes a random variable that is equal in distribution to J 1 . Our first objective is to find the distribution of the number of jobs in the system at time t, in the sequel denoted by N (t). This can be found in several ways; because of the appealing underlying intuition, we here provide an argument in which we approximate the arrival rate on intervals of length by a constant, and then let ↓ 0. This procedure works as follows. We let (t) = (ω, t) be an arbitrary sample path of the driving shot-noise process. Given (t), the number of jobs that arrived in the interval [k , (k + 1) ) and are still in the system at time t has a Poisson distribution . . are i.i.d. standard uniform random variables. Summing over k yields that the number of jobs in the system at time t has a Poisson distribution with parameter The argument above is not new: a similar observation was mentioned in, for example [6], for deterministic rate functions. Since (·) is actually a stochastic process, we conclude that the number of jobs has a mixed Poisson distribution, with the expression in Eq. (3) as random parameter. As a consequence, we find by conditioning on F t , We have found the following result.

Theorem 3.1 Let (·) be a shot-noise process. Then
Proof The result follows directly from Lemma 2.1 and Eq. (4).
In Theorem 3.1 we found that N (t) has a Poisson distribution with the random parameter given in Eq. (3). This leads to the following expression for the expected value: In addition, by the law of total variance we find The latter expression we can further evaluate, using an approximation argument that resembles the one we used above. Using a Riemann sum approximation, we find . We thus have that (7) equals We can make this more explicit using the corresponding formulas in (2). Example 3.2 (Exponential case) Consider the case in which J is exponentially distributed with mean 1/μ and (0) = 0. Then we can calculate the mean and variance explicitly. For μ = r , where the function h r,μ (·) is defined by For the variance, we thus find for μ = r and for μ = r

Asymptotic analysis
This subsection focuses on deriving a functional central limit theorem (FCLT) for the model under study, after having appropriately scaled the shot rate of the shot-noise process. In the following, we assume that the service requirements are exponentially distributed with rate μ, and we point out how it can be generalized to a general distribution in Remark 3.6 below. We follow the standard approach to derive the FCLT for infinite-server queueing systems; we mimic the argumentation used in, for example [1,14]. As the proof has a relatively large number of standard elements, we restrict ourselves to the most important steps.
We apply a linear scaling to the shot rate of the shot-noise process, i.e., ν → nν. It is readily checked that under this scaling, the steady-state level of the shot-noise process and the steady-state number of jobs in the queue blow up by a factor n. It is our objective to prove that, after appropriate centering and normalization, the process recording the number of jobs in the system converges to a Gaussian process.
In the n-th scaled model, the number of jobs in the system at time t, denoted by N n (t), has the following (obvious) representation: with A n (t) denoting the number of arrivals in [0, t], and D n (t) the number of departures, Here, A n (t) corresponds to a Cox process with a shot-noise-driven rate, and therefore we have, with n (s) the shot-noise in the scaled model at time s and S A (·) a unit-rate Poisson process, in line with our previous assumptions, we put n (0) = 0. For our infinite-server model the departures D n (t) can be written as, with S D (·) a unit-rate Poisson process (independent of S A (·)), We start by identifying the average behavior of the process N n (t). Following the reasoning of [1], assuming that N n (0)/n ⇒ ρ(0) (where '⇒' denotes weak convergence), N n (t)/n converges almost surely to the solution of Now we move to the derivation of the FCLT. Following the approach used in [1], we proceed by studying an FCLT for the input rate process. To this end, we first definê The following lemma states thatK n (·) converges to an integrated Ornstein-Uhlenbeck (OU) process, corresponding to an OU processˆ (·) with a speed of mean reversion equal to r , long-run equilibrium level 0, and variance σ 2 := ν E B 2 /(2r ).

Lemma 3.3 Assume that for the shot sizes, distributed as B, it holds that
in whichˆ satisfies, with W 1 (·) a standard Brownian motion, Proof This proof is standard; for instance from [2,Prop. 3], by putting the λ d in that paper to zero, it follows thatˆ n (·) ⇒ˆ (·). This impliesK n (·) ⇒K (·), using (10) together with the continuous mapping theorem.
Interestingly, the above result entails that the arrival rate process displays meanreverting behavior. This also holds for the job count process in standard infinite-server queues. In other words, the job count process in the queueing system we are studying can be considered as the composition of two mean-reverting processes. We make this more precise in the following.
From now on we consider the following centered and normalized version of the number of jobs in the system: We assume thatN n (0) ⇒N (0) as n → ∞. To prove the FCLT, we rewriteN n (t) in a convenient form. Mimicking the steps performed in [1] or [14], withS and using the relation (9), we eventually obtain Our next goal is to apply the martingale FCLT to the martingales R n (t)/ √ n; see, for background on the martingale FCLT, for instance [7] and [17]. The quadratic variation equals Appealing to the martingale FCLT, the following FCLT is obtained.
with W 2 (·) a standard Brownian motion that is independent of the Brownian motion W 1 (·) we introduced in the definition ofK (·).
Remark 3.5 In passing, we have proven that the arrival process as such obeys an FCLT. WithÂ we find thatÂ n (t) ⇒Â(t) as n → ∞, wherê the last equality follows from the fact that ρ(·) satisfies (9).

Remark 3.6
The FCLT can be extended to non-exponential service requirements, by making use of [15,Thm. 3.2]. Their approach relies on two assumptions: • The arrival process should satisfy an FCLT; • The service times are i.i.d. nonnegative random variables with a general c.d.f. independent of the arrival process.
As noted in Remark 3.5, the first assumption is satisfied for the model in this paper. The second assumption holds as well. In the non-exponential case, the results are less clean; in general, the limiting process can be expressed in terms of a Kiefer process, cf. for example [4].

Networks
Now that the reader is familiar with the one-dimensional setting, we extend this to networks. In this section, we focus on feedforward networks in which each node corresponds to an infinite-server queue. Feedforward networks are defined as follows. E) be a directed graph with nodes V and edges E. The nodes represent infinite-server queues and the directed edges between the facilities demonstrate how jobs move through the system. We suppose that there are no cycles in G, i.e., there is no sequence of nodes, starting and ending at the same node, with each two consecutive nodes adjacent to each other in the graph, consistent with the orientation of the edges.
We focus on feedforward networks to keep the notation manageable. In Theorem 4.6, we derive the transform of the numbers of jobs in all nodes, jointly with the shot-noise process(es) for feedforward networks. Nonetheless, we provide Example 4.5, to show that analysis is in fact possible if there is a loop, but at the expense of more involved calculations.
Since all nodes represent infinite-server queues, one can see that whenever a node has multiple input streams, it is equivalent to multiple infinite-server queues that work independently from each other, but have the same service speed and induce the same service requirement for arriving jobs. Consider Fig. 3 for an illustration. The reason why this holds is that different job streams move independently through the system, without creating waiting times for others. Therefore, merging streams do not increase the complexity of our network. The same holds for 'splits' in job streams. By this we mean that after jobs finished their service at a server, they move to server i with probability q i (with i q i = 1). Then, one can simply sample the entire path that the job will take through the system, at the arrival instance at its first server.
If one recognizes the above, then all feedforward networks reduce to parallel tandem systems in which the first node in each tandem system is fed by external input. The procedure to decompose a network into parallel tandems consists of finding all paths between nodes in which jobs either enter or leave the system. Each of these paths will subsequently be considered as a tandem queue, which are then set in parallel. To build up to the main result, we first study tandem systems in Sect. 4.1. Subsequently, we put the tandem systems in parallel in Sect. 4.2 and finally we present the main theorem and some implications in Sect. 4.3.

Tandem systems
As announced, we proceed by studying tandem systems. In Sect. 4.2 below, we study d parallel tandem systems. In this subsection, we consider the i -th of these tandem systems, where i = 1, . . . , d.
Suppose that tandem i has S i service facilities and the input process at the first node is Poisson, with a shot-noise arrival rate i (·). We assume that jobs enter node i1. When they finish service, they enter node i2, etc., until they enter node i S i after which they leave the system. We use i j as a subscript referring to node j in tandem system i and we refer to the node as node i j. Hence N i j (t) and J i j denote the number of jobs in node i j at time t, and a copy of a service requirement, respectively, where j = 1, . . . , S i .
Fix some time t > 0. Again we derive results by splitting time into intervals of length . Denote by M i j (k, ) the number of jobs present in node i j at time t that have entered node i1 between time k and (k + 1) ; as we keep t fixed we suppress it in our notation. Because jobs are not interfering with each other in the infinite-server realm, we can decompose the transform of interest: Supposing that the arrival rate is a deterministic function of time λ i (·), by conditioning on the number of arrivals in the k-th interval, in which p i (u) ( p i j (u), respectively) denotes the probability that the job that entered tandem i at time u has already left the tandem (is in node j, respectively) at time t.
Note that Recognizing a Riemann sum and letting ↓ 0, we conclude that Eq. (12) takes the following form: In case of a stochastic rate process i (·), we obtain Therefore, it holds that and we consequently find

Parallel (tandem) systems
Now that the tandem case has been analyzed, the next step is to put the tandem systems as described in Sect. 4.1 in parallel. We assume that there are d parallel tandems. There are different ways in which dependence between the parallel systems can be created. Two relevant models are listed below, and illustrated in Fig. 4.  In Model (M1), correlation between the shot-noise arrival rates induces correlation between the numbers of jobs in the different queues. In Model (M2), correlation clearly appears because all tandem systems have the same input process. Of course, the tandem systems will not behave identically because the jobs may have different service requirements. In short, correlation across different tandems in Model (M1) is due to linked arrival rates, and correlation in Model (M2) is due to simultaneous arrival epochs. We feel that both versions are relevant, depending on the application, and hence we analyze both.
Analysis of (M1)-Suppose that the dependency is of the type as in Model (M1). This means that the shots, in each component of , occur simultaneously. Recall the definition of f i as stated in Eq. (13). It holds that where the last equality holds due to (15).
Analysis of (M2)-Now suppose that the dependency in this model is of type (M2), i.e., there is one shot-noise process that generates simultaneous arrivals in the parallel tandem systems. First we assume a deterministic arrival rate function λ(·). Let M i j (k, ) be the number of jobs present in tandem system i at node j at time t that have arrived in the system between k and (k + 1) . Note that To further evaluate the right-hand side of the previous display, we observe that we can write in this definition p 1 ,..., d ≡ p 1 ,..., d (u) equals the probability that a job that arrived at time u in tandem i is in node i at time t [cf. Eq. (14)]. The situation that i = S i + 1 means that the job left the tandem system; we define z i,S i +1 = 1.
Remark 4.4 (Routing) Consider a feedforward network with routing. As argued in the beginning of this section, the network can be decomposed as a parallel tandem system. In case there is splitting at some point, then one decomposes the network as a parallel system, in which each tandem i receives the job with probability q i , such that q i = 1. This can be incorporated simply by adjusting the probabilities contained in f i in Eq. (16), which are given in Eq. (14), so that they include the event that the job joined the tandem under consideration. For instance, the expression for p i (u) in the left equation in (14) would become where Q is a random variable with a generalized Bernoulli (also called 'categorical') distribution, where P(job is assigned to tandem i) = P(Q = i) = q i , for i = 1, . . . d, with q i = 1; the right equation in (14) is adjusted similarly. Other than that the analysis is the same for the case of splits.
Remark 4.5 (Networks with loops) So far we have only considered feedforward networks. Networks with loops can be analyzed as well, but the notation becomes quite cumbersome. To show the method in by which networks with loops and routing can be analyzed, we consider a specific example. Suppose that arrivals enter node one, after which they enter node two. After they have been served in node two, they go back to node one with probability η, or leave the system with probability 1 − η. In this case, with similar techniques as before, we can find in which job(u) is the job that arrived at time u and we are examining the system at time t. Now, if we denote service times in the j-th node by J ( j) , then, at a specific time t, Analogously, P( job(u)is in node 1) equals, by conditioning on the job having taken k loops, For example, in the case where all J ( j) i are independent and exponentially distributed with mean 1/μ, we can calculate those probabilities explicitly. Indeed, if we denote by Y a Poisson process with rate μ, then, for example A similar calculation can be done for the probability that the job is in node one.
Recalling that a sum of exponentials has a Gamma distribution, we can write where F (2m+2,μ) denotes the distribution function of a -distributed random variable with rate μ and shape parameter 2m + 2.

Main result
In this subsection, we summarize and conclude with the following main result.

and where g(s, v) is a vector-valued function in which component i is given by
with f (·, ·) as defined in Eq. (17).
Next we calculate covariances between nodes in tandem and parallel thereafter.
Covariance in Tandem System-Consider a tandem system consisting of two nodes and we want to analyze the covariance between the number of jobs in the nodes. Dropping the index of the tandem system, denote by N 1 (·) and N 2 (·) the number of jobs in nodes 1 and 2, respectively. Using Eq. (16), differentiation yields cf. Eq. (2) for the last equality.
Covariance parallel (M1)-Consider a parallel system consisting of two nodes only.
In order to study covariance in the parallel (M1) case, we need a result about the covariance of the corresponding shot-noise process.
Lemma 4.7 Let 1 (·), 2 (·) be shot-noise processes of which the jumps occur simultaneously according to a Poisson arrival process with rate ν. Let the decay be exponential with rate r 1 , r 2 , respectively. Then it holds that, for δ > 0, which, in the case 1 = 2 , reduces to corresponding to [16, p. 394] and Eq. (2).
Proof See Appendix 2.
By making use of Eq. (16), we find

This implies
where we made use of the fact that, for u ≥ v, cf. Lemma 4.7.
Covariance parallel (M2)-Extracting the mixed moment from the transform in Eq. (18), we derive directly that

This implies
The following proposition compares the correlations present in Model (M1) and (M2).
In the proposition, we refer to the number of jobs in queue j in Model (Mi) at time t as N (i) j (t), for i = 1, 2. We find the anticipated result that the correlation in Model (M2) is stronger than in Model (M1).
Proof Because of the assumption 1 The expressions for the covariances, which are derived earlier in this section, imply that Cov N which is nonnegative, as desired.

Concluding remarks
We have considered networks of infinite-server queues with shot-noise-driven Coxian input processes. For the single queue, we found explicit expressions for the Laplace transform of the joint distribution of the number of jobs and the driving shot-noise arrival rate, as well as a functional central limit theorem of the number of jobs in the system under a particular scaling. The results were then extended to a network context: we derived an expression for the joint transform of the numbers of jobs in the individual queues, jointly with the values of the driving shot-noise processes.
We included the functional central limit theorem for the single queue, but it is anticipated that a similar setup carries over to the network context, albeit at the expense of considerably more involved notation. Our future research will include the study of the departure process of a single queue; the output stream should remain Coxian, but of another type than the input process. recognizing a Riemann sum, With P B (t) a Poisson process with rate ν and the U i i.i.d. samples from a uniform distribution on [0, 1], it holds that We thus obtain that the expression in (20) equals (where the equality follows by interchanging the order of the summations) which behaves as Furthermore, we have the representation We conclude that E z N (t) e −s (t) equals Conditioning on the values of P B ( ) − P B (( −1) ), for = 1, . . . , t/ , and using that the B i are i.i.d., we find that the expression in Eq. (21)) is equal to the limit, as ↓ 0, of the following expression The lemma now follows from continuity of the exponent and the definition of the Riemann integral.

Appendix 2: Proof of Lemma 4.7
Let P B (·) be the Poisson process with rate ν, corresponding to the occurrences of shots, and let E t,δ (n) be the event that P B (t + δ) − P B (t) = n. By conditioning on the number of shots in the interval (t, t + δ], we find We proceed by rewriting the conditional expectation as ..,t n ,δ (n)) dt 1 . . . dt n , denoting by F t 1 ,...,t n ,δ (n) the event E t,δ (n) and the arrival epochs are t 1 , . . . , t n . Note that we have, due to Eq. (1), conditional on F t 1 ,...,t n ,δ (n), the distributional equality 2 (t + δ) = 2 (t)e −r 2 δ + n i=1 B i2 e −r 2 (t+δ−t i ) , and consequently E 1 (t) 2 (t + δ) | F t 1 ,...,t n ,δ (n) Note that for all i = 1, . . . , n we have After unconditioning Eq. (23) with respect to the arrival epochs by integrating over all t i from t to t + δ and dividing by δ n , we thus obtain (1 − e −r 2 δ )n E B 12 and hence, denoting i := lim t→∞ i (t) for i = 1, 2, where we made use of E i = ν E B 1i /r i . It follows that Cov( 1 (t), 2 (t + δ)) = E 1 (t) 2 where the equality E 2 (t + δ) = E 2 (t)e −r 2 δ + (1 − e −r 2 δ ) E 2 is used, which can be directly checked using the expressions for the mean in Eq. (2). This proves the first equality in Eq. (19). The proof of the second equality follows from