Efficient computation of overlapping variance estimators for simulation.

For a steady-state simulation output process, we formulate efficient algorithms to compute certain estimators of the process variance parameter (i.e., the sum of covariances at all lags), where the estimators are derived in principle from overlapping batches separately and then averaged over all such batches. The algorithms require order-of-sample-size work to evaluate overlapping versions of the area and Cramer-von Mises estimators arising in the method of standardized time series. Recently, Alexopoulos et al. showed that, compared with estimators based on nonoverlapping batches, the estimators based on overlapping batches achieve reduced variance while maintaining similar bias asymptotically as the batch size increases. We provide illustrative analytical and Monte Carlo results for M/M/1 queue waiting times and for a first-order autoregressive process. We also present evidence that the asymptotic distribution of each overlapping variance estimator can be closely approximated using an appropriately rescaled chi-squared random variable with matching mean and variance.

Key words: steady-state simulation; simulation output analysis methods; method of batch means; method of standardized time series; area variance estimator; Cramer-von Mises variance estimator; nonoverlapping variance estimator; overlapping variance estimator

1. Introduction

The simplest goal in a steady-state simulation experiment is often to estimate the unknown mean [mu] of a resulting output process, {[Y.sub.i]: i = 1,2,...}, which is taken to be stationary. Based on a simulation-generated time series of length n, the obvious estimator for [mu] is the sample mean, [bar.Y.sub.n] = [n.sup.-1] [[summation].sub.i=1.sup.n] [Y.sub.i]. In addition, the careful experimenter ought to provide an estimate for the variance of the sample mean, Var([bar.Y.sub.n]) = E[([bar.Y.sub.n] - [mu])[.sup.2]]. To characterize the asymptotic behavior of [bar.Y.sub.n] as the run length (sample size) n tends to infinity, the experimenter could estimate either (i) the quantities {[[sigma].sub.n.sup.2] [equivalent to] nVar ([bar.Y.sub.n]): n = 1,2,...}; or (ii) the variance parameter, [[sigma].sup.2] [equivalent to] [lim.sub.n[right arrow][infinity]] [[sigma].sub.n.sup.2], provided that this limit exists. Among the many techniques in the literature for analyzing steady-state simulation output and for estimating [[sigma].sup.2] are the following: nonoverlapping batch means (NBM) (Schmeiser 1982); overlapping batch means (OBM) (Meketon and Schmeiser 1984); and standardized time series (STS) (Schruben 1983). The STS technique gives rise to such estimators as the STS area estimator (Goldsman et al. 1990) and the STS Cramer-von Mises (CvM) estimator (Goldsman et al. 1999).

A common strategy used with the above estimators involves batching the observations. For instance, the NBM, STS area, and STS CvM estimators split the set of observations into adjacent disjoint (nonoverlapping) batches, where each batch is a subseries (i.e., a group of consecutive observations) that is presumably large enough to provide a representative "snapshot" of the overall dependency structure of the entire simulation-generated time series {[Y.sub.i]: i = 1,..., n}. In the case of NBM, we then assume that the resulting sample (batch) means computed from each batch are approximately independent and identically distributed (i.i.d.) normal random variables; and finally we appropriately rescale the sample variance of the batch means so as to obtain an estimator of [[sigma].sup.2]. In the case of the STS estimators, we compute a separate STS estimator of [[sigma].sup.2] from each batch; we assume the resulting STS estimators are i.i.d.; and then we average those estimators to obtain an overall batched estimator of [[sigma].sup.2].

On the other hand, the OBM method uses overlapping batches, so that the corresponding batch means are not independent (though they are identically distributed and asymptotically normal). As in the computation of the NBM variance estimator, the OBM estimator of [[sigma].sup.2] is an appropriately rescaled version of the sample variance of the overlapping batch means; and it can be shown that as the batch size m [right arrow] [infinity] (so that the simulation run length n [right arrow] [infinity]), the OBM estimator asymptotically has (nearly) the same bias as, but smaller variance than, the corresponding NBM estimator for [[sigma].sup.2]. Alexopoulos et al. (2006b) formulated overlapping versions of the STS area and CvM estimators; and they found that as the batch size m [right arrow] [infinity], the STS overlapping estimators asymptotically had (nearly) the same bias as, but substantially smaller variance than, their counterparts computed from nonoverlapping batches. In this paper, we present numerically efficient algorithms for computing the new variance estimators, and we evaluate the statistical performance of these estimators in several analytical and Monte Carlo examples.

The rest of this paper is organized as follows. First, in Section 2 we provide a self-contained synopsis of relevant background material. In Section 3 we formulate numerically efficient algorithms for computing the proposed overlapping variance estimators, and we show these algorithms to possess the unexpected advantage of requiring execution times that are linear in the sample size. To complement the theoretical results presented in Alexopoulos et al. (2006b) on the large-sample moment and distributional properties of our new variance estimators, in Section 4 we present both theoretically exact and empirically estimated small-sample results for some illustrative examples with finite but progressively increasing batch sizes. Moreover, Section 4 provides both theoretical and empirical evidence that for a sufficiently large batch size, the distribution of each overlapping variance estimator can be closely approximated by the distribution of a rescaled chi-squared random variable whose scaling factor and degrees of freedom are set to match the mean and variance of the target distribution. To round out the development in Section 4, we formulate approximate confidence intervals for [mu] and [[sigma].sup.2] based on our approximate distributions for the overlapping variance estimators. Finally, in Section 5 we summarize the main contributions of this work and make recommendations for future research.