{'fraction_non_alphanumeric': 0.06043802162101607, 'fraction_numerical': 0.012714949868369461, 'mean_word_length': 4.148212226066898, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Measuring the Initial Transient: Reflected Brownian Motion

We analyze the convergence to equilibrium of one-dimensional reflected Brownian motion (RBM) and compute a number of related initial transient formulae. These formulae are of interest as approximations to the initial transient for queueing systems in heavy traffic, and help us to identify settings in which initialization bias is significant. We conclude with a discussion of mean square error for RBM. Our analysis supports the view that initial transient effects for RBM and related models are typically of modest size relative to the intrinsic stochastic variability, unless one chooses an especially poor initialization.


Introduction
This paper is concerned with using one-dimensional reflected Brownian motion (RBM) as a theoretical vehicle for studying the initial transient problem. Given that RBM is a commonly used approximation to a wide variety of different queueing models, the initial transient behavior of RBM can be viewed as being representative of a large class of simulation models in which congestion is a key factor.
The stochastic process X = (X(t) : t ≥ 0) is said to be a (one-dimensional) RBM if it satisfies the stochastic differential equation (SDE) where B = (B(t) : t ≥ 0) is standard Brownian motion and L is a continuous nondecreasing process that increases only when X is at the origin (so that I(X(t) > 0)dL(t) = 0). In particular, the process L is a "boundary process" that serves to keep X nonnegative as befits an approximation to a queue. The parameter −r represents the "drift" of the RBM, and σ is its "volatility" parameter.
To illustrate the sense in which RBM can be used to approximate a queue, consider a system with a single queue that is being fed by a renewal arrival process in which χ denotes a generic interarrival time random variable (rv). Customers are served by one of m identical servers, in the order in which they arrive. The service times are independent and identically distributed (iid) across the servers and across the customers, and are also independent of the interarrival times. If V is a rv having the common service time distribution, set λ = 1/Eχ, µ = 1/EV and σ 2 A = var χ, σ 2 S = var V . It is well known that if Z(t) is the numberin-system at time t, then where X is an RBM with r = mµ − λ and σ 2 = λ 3 σ 2 A + mµ 3 σ 2 S , provided that mµ − λ is small (so that the system is in "heavy traffic") and t is of the order (mµ − λ) −2 . Here, D ≈ means "has approximately the same distribution as", and the rigorous support rests on a so-called "heavy traffic" limit theorem; see, for example, Iglehart and Whitt (1970).
This paper begins with a discussion of the convergence to equilibrium of RBM (Section 2). It subsequently develops closed-form expressions for various initial transient quantities associated with RBM, distinguishing between "functional" (Section 3) and "distributional" (Section 4) perspectives. These expressions can be used to help plan steady-state/equilibrium simulations of queueing models that can be approximated by RBM as well as to identify settings in which initial bias is significant (Section 5). We conclude by deriving a decomposition of mean square error (MSE) in the setting of RBM (Section 6).
Though only in the context of RBM, the results in this paper are intended develop general insights into the initial transient problem that can be of potentially broader applicability.

Convergence to Equilibrium for RBM
It is well known that if r > 0, then X has an equilibrium, in the sense that σ 2 (see, for example, Harrison (1985) pg.94). A key question in the study of the initial transient problem for X is its rate of convergence to equilibrium. One vehicle for studying this question is the well-known formula (2.1) for the transition density of X; see Harrison (1985) pg.49. Here, p(t, x, y)dy ∆ = P (X(t) ∈ dy | X(0) = x), Φ(x) = P (N (0, 1) ≤ x) (where N (0, 1) denotes a normal rv with mean 0 and unit variance), and φ(x) is the density associated with φ. But it is difficult to "read off" the rate of convergence from (2.1).
However, an alternative representation for the transition density of X can be computed. Recall that the rate at which the transition probabilities of a Markov jump process converge to their respective equilibrium probabilities can easily be determined once the eigenvalues and eigenvectors of the rate matrix are known. Something similar can be implemented in the RBM setting. This leads to an alternative representation of the transition density known as the spectral decomposition.
To determine the eigenvalues and eigenfunctions associated with X, note that Itô's formula guarantees that d(e −λt u(X(t))) = −λe −λt u(X(t))dt + e −λt u (X(t))dX(t) + e −λt u (X(t)) 2 where L is the second order differential operator given by Note that because L increases only when X is at the origin, u (X(t))dL(t) = u (0)dL(t). Consequently, if u (0) = 0 and the stochastic integral is integrable, e −λt u(X(t)) − t 0 e −λs ((L u)(X(s)) − λu(X(s)))ds is a martingale. It follows that if u also satisfies L u = λu, then E x u(X(t)) = e λt u(x) for t, x ≥ 0, where E x (·) ∆ = E(· | X(0) = x). In other words, u is an eigenfunction of L (and X), and λ is its associated eigenvalue. Since the eigenvalues of a Markov process automatically have non-positive real parts, sending t → ∞ allows us to conclude that Eu(X(∞)) = 0, and the rate of convergence of E x u(X(t)) is exponentially fast with associated rate parameter λ.
Remark: This makes clear that the rate of convergence of E x f (X(t)) to Ef (X(∞)) can depend on the choice of f . In particular, when f is an eigenfunction, the rate depends on its associated eigenvalue.
Eigenfunction/eigenvalue pairs for RBM exist only when λ = 0 and λ ≤ −γ, where γ ∆ = r 2 /(2σ 2 ). Any λ ∈ (−∞, −γ] is an eigenvalue, with associated eigenfunction where s(λ) = σ −2(λ + γ). As a consequence, the distance between the top two eigenvalues of L , namely 0 and −γ, is equal to γ. This distance is known, in the Markov process literature, as the spectral gap of X. Given the eigenvalues and eigenfunctions, there is a standard "recipe" for constructing the transition density for reversible diffusion processes (of which RBM is one) that can be found, for example, on pg.332 of Karlin and Taylor (1981). When specialized to the RBM setting, one obtains the spectral decomposition for p(t, x, y), namely (2005)). We thus find that, for f appropriately integrable, The spectral representation (2.2) makes clear that the spectral gap γ is precisely the exponential rate constant governing the rate at which E x f (X(t)) converges to Ef (X(∞)). Furthermore, when t is large, it is primarily the "projection" of f onto u −γ (i.e. the magnitude of f, u λ ) that determines the magnitude of E x f (X(t))−Ef (X(∞)) (given that the integral is largely determined by the integrand's contribution from a neighborhood of λ = −γ).

The Initial Transient Effect
Given a performance measure f , we have studied in Section 2 the rate of convergence of E x f (X(t)) to its equilibrium value Ef (X(∞)). Our goal here is to compute the magnitude of the initial transient effect, assuming that no deletion is implemented. This can inform our decision as to how serious the initial transient effect is, and whether/how deletion is warranted. Because the typical estimator used to compute an equilibrium quantity in a simulation context is a time-average, the effect of the initial transient in the steady-state simulation setting is, in some sense, an integrated version of the theory of Section 2.
To compute the effect of the initial transient, note that if h is twice continuously differentiable with h (0) = 0, then Itô's formula yields

It follows that if the stochastic integral is integrable and
In view of our discussion in Section 2, we expect that Hence, the constant h c (x) expresses (up to an exponentially small order) the bias of the estimator t −1 t 0 f (X(s))ds induced by the initial transient.
We now turn to computing h c (x). The solution h to Poisson's equation Given the performance measure f , the functional unconditional initial transient effect (UNITE) measure is given by for a given initial distribution µ, while the functional conditional initial transient effect (CITE) measure is defined byβ where δ x (·) is a unit point mass distribution at x; see Wang and Glynn (2014) for a more detailed discussion of these measures. Note that in a single replication setting, the initial transient effect is determined by the random placement of X(0), so that averaging the effect over the initial distribution µ (as in CITE) seems reasonable. Given this viewpoint, it is then appropriate to viewβ f (π) as a benchmark against whichβ f (µ) for other initializations µ can be compared. (After all, initializing with π makes X a stationary process in which no initial transient is present.) In particular, we can now separate the state space into "good states" and "bad states", depending on whether β f (δ x ) ≤ cβ f (π) for the state x or not, where c is a given constant. For example, for c = 1, the good states are all those states x for which β f (δ x ) is smaller than the CITE measure associated with π. Of course, the set of good states is sensitive to the choice of c. Hence, it is of interest to study the dependence of the set of good states on the parameter c; see Figure 1 below. We note that the notion of good state is related to the "typical state" idea introduced by Grassmann (2011) and Grassmann (2014).
For f (x) = x, the set of good states G (r, σ 2 ; c) For this performance measure, G (r, σ 2 ; c) = EX(∞)G (1, 2; c), so it suffices to graph G (r, σ 2 ) ∆ = {(x, c) : x ∈ G (r, σ 2 ; c), c ≥ 0} only for EX(∞) = 1; the resulting graph can be found below. Note that for c = 1, the set of good states already covers [0, 2.09], so that the set of initial values providing single-replication bias characteristics roughly comparable to that associated with initializing under π is quite robust (in fact, more than twice the mean).

The Distributional Initial Transient
The convergence of X(t) to X(∞) involves the entire distribution of X (as opposed to only the convergence of E x f (X(t)) for a single performance measure f ), so that much of the probability literature is concerned with the rate of convergence at which the distribution of X(t) approaches that of X(∞). In fact, one commonly used measure for judging the rate of convergence to equilibrium is the weighted total variation norm defined by sup |f |≤w |E µ f (X(t)) − Ef (X(∞))| for a given weight function w : R + → R + .
By analogy, with the discussion of Section 3, it is natural to study the distributional UNITE measure given by and the distributional CITE measure defined bỹ In view of the calculation of Section 3, specifically (3.1), it is evident that According to Pollard (2002) pg.60, Hence, if w(y) = y p for p ≥ 0, So, as in Section 3, we can analogously define the set of "good states" as so that H (r, σ 2 ; c) again scales similarly as does the set G (r, σ 2 ; c) of Section 3, due to the homogeneity of w(y) = y p . (For the function w(y) = e θy , such a scaling relationship does not hold.) By numerically computingβ(π) (via a numerical integration of β(δ x ) against π), the graph of the set H (1, 2) H (1, 2; c), c ≥ 0} can be determined. In particular, it is computed below for w(x) ≡ 1; see Figure 2. Unlike the set G (1, 2) of Section 3, it does not touch the x-axis. This is not surprising, because in the functional setting of Section 3, it will often be the case that the bias is monotone in the initial state. For example, for nondecreasing f , the stochastic monotonicity of RBM ensures that the bias will be negative for small values of x, and positive for large values of x, yielding (using a continuity argument) the existence of an intermediate x at which β f (δ x ) = 0. On the other hand, in the distributional setting, β(δ x ) involves looking at the largest possible bias over a large class of functions f , and there typically will be no initial x at which β(δ x ) will vanish.
We further note that the shape of H (1, 2) is quite different than that of G (1, 2). For example, the minimizer of β(δ x ) as a function of x is achieved at a point x * that is smaller than that of β f (δ x ) for f (x) = x. This occurs because β(δ x ) is defined in terms of the bias of bounded functions, while the set G (1, 2) of Section 3 was computed for the identity mapping (which is unbounded). In computing the expectation of such a performance measure, it is advantageous to initialize the process at a larger value, because such an initialization will lead to higher likelihood paths that will quickly sample the large state values that typically contribute the most to the expectation of the performance measure. We further note that the set of good x's associated with H (1, 2) is a bit smaller than that of Section 3 (in part, because H (1, 2) involves a worst case bias, where as G (1, 2) is determined only by a single function f ). Nevertheless, the set of good states, with values of β(δ x ) of a magnitude less than or equal to that associated with initializing under the equilibrium distribution (i.e., β(π)), is large, and includes all the states x with a value less than or equal to 1.84 EX(∞).
On the other hand, the systematic error in α(t) (namely, the bias) was analyzed in Section 3 and determined to be β f (δ x )/t for t large, assuming that the process was initialized at state x. It is clear that for t large enough, the stochastic variability dominates the bias. Specifically, for t ≥ t * (x) ∆ = β f (δ x ) 2 π/(2κ 2 ), the bias contribution is smaller than the error due to stochastic variability.
To help put the quantity t * (x) in perspective, note that the CLT suggests that a relative error * = 2κ(2/π) 1/2 /(t * (x) 1/2 |Ef (X(∞))|) is achieved at a run-length t * (x) (contributed in equal measures from stochastic variability and bias). In other words, * (x) is the relative precision at which a simulation designed to achieve such an error tolerance will have an error that is contributed equally from stochastic variability and initial transient bias. At all smaller values of the relative precision, the stochastic variability dominates the error. In particular, to achieve a relative error tolerance of = ν * (x) (with ν < 1), the associated run-length t that must be used is such that the error due to initial transient bias at such a run-length is roughly a proportion ν of the stochastic error. Thus, the quantity * (x) is an important measure of the threshold error tolerance at which the initial transient is no longer a dominant source of error in computing steady-state quantities. Figure 3 below provides a graph of the threshold error tolerance as a function of x for f (y) = y, with (r, σ 2 ) = (1, 2). Given that a relative error precision of 10% or less is typically desired, we note that the set of states x for which * (x) ≥ 0.1 is large, specifically the interval [0, 10.19].
is a deterministic function tending to 0 as t → ∞, and (X * , B) = ((X * (t), B(t)) : −∞ < t < ∞) is such that X * is a stationary RBM driven by the Brownian motion B. Of course, Since X * is a reversible one-dimensional diffusion driven by B (see, for example, Kent (1978)), where D = denotes equality in distribution. Consequently, (6.1) equals because of the martingale property of the stochastic integral. Furthermore, E x h 2 c (X(t)) = Eh 2 c (X(∞))+o(1) as t → ∞ and E x h c (X(t)) = Eh c (X(∞))+ o(1) as t → ∞. Finally, by following the same argument as in Section 3, we find that as t → ∞, where k c is the solution of the Poisson's equation associated with h 2 c , namely k c satisfies (L k c )(x) = −(h c (x) 2 − Eh c (X(∞)) 2 ), x ≥ 0 subject to k c (0) = 0, Ek c (X(∞)) = 0. In view of (3.1), With k c now computed, the mean square error of α(t) is given by as t → ∞. In particular, for f (x) = x, σ 6 2r 4 t + 6r 4 x 4 + 8r 3 σ 2 x 3 + 6r 2 σ 4 x 2 − 3σ 8 24r 6 t 2 + 1 t 2 o (1) as t → ∞. The term h c (x) 2 /t 2 is the squared bias contribution to the MSE due to the initial transient. The expression (6.2) makes clear that the MSE includes other state-dependent contributions of the same order of magnitude (that are contributed by the variance of α(t) rather than the bias), namely σ 2 k c (x)/t 2 . Hence, a full analysis of the MSE impact of the effect of the initial transient should also (ideally) include an analysis of this variance term involving k c , in addition to the bias contribution that is typically included in such an MSE analysis.

Conclusion
We have developed various formulae related to the initial transient problem for RBM. These formulae can be used directly in a simulation context, to approximate the impact of the initial transient for queueing simulations for which RBM is a suitable guide (eg. simulations of queues in heavy traffic). Our formulae also make clear a key insight that is likely true in a much broader class of simulations. In particular, for RBM, there is typically a robust set of initializing states for which the impact of the initial transient is roughly comparable to that associated with initializing in equilibrium. Since initializing in equilibrium corresponds to a setting in which there is no initial transient, this suggests that one should not worry excessively about the initial transient unless one has inadvertently initialized the simulation with a very poor ("bad") choice of state. If one instead initializes with a reasonable ("good") choice of state, the key element to a successful calculation of Ef (X(∞)) is ensuring that the run-length is long enough to ensure that the stochastic variability has been reduced to a level commensurate with the desired accuracy. This suggests that a focus of initial transient