Input uncertainty is an aspect of simulation model risk that arises when the driving input distributions are derived or “fit” to real-world, historical data. Although there has been significant progress on quantifying and hedging against input uncertainty, there has been no direct attempt to reduce it via better input modeling. The meaning of “better” depends on the context and the objective: Our context is when (a) there are one or more families of parametric distributions that are plausible choices; (b) the real-world historical data are not expected to perfectly conform to any of them; and (c) our primary goal is to obtain higher-fidelity simulation output rather than to discover the “true” distribution. In this paper, we show that frequentist model averaging can be an effective way to create input models that better represent the true, unknown input distribution, thereby reducing model risk. Input model averaging builds from standard input modeling practice, is not computationally burdensome, requires no change in how the simulation is executed nor any follow-up experiments, and is available on the Comprehensive R Archive Network (CRAN). We provide theoretical and empirical support for our approach.

Inference-based optimization via simulation, which substitutes Gaussian process (GP) learning for the structural properties exploited in mathematical programming, is a powerful paradigm that has been shown to be remarkably effective in problems of modest feasible-region size and decision-variable dimension. The limitation to “modest” problems is a result of the computational overhead and numerical challenges encountered in computing the GP conditional (posterior) distribution on each iteration. In this paper, we substantially expand the size of discrete-decision-variable optimization-via-simulation problems that can be attacked in this way by exploiting a particular GP—discrete Gaussian Markov random fields—and carefully tailored computational methods. The result is the rapid Gaussian Markov Improvement Algorithm (rGMIA), an algorithm that delivers both a global convergence guarantee and finite-sample optimality-gap inference for significantly larger problems. Between infrequent evaluations of the global conditional distribution, rGMIA applies the full power of GP learning to rapidly search smaller sets of promising feasible solutions that need not be spatially close. We carefully document the computational savings via complexity analysis and an extensive empirical study. Summary of Contribution: The broad topic of the paper is optimization via simulation, which means optimizing some performance measure of a system that may only be estimated by executing a stochastic, discrete-event simulation. Stochastic simulation is a core topic and method of operations research. The focus of this paper is on significantly speeding-up the computations underlying an existing method that is based on Gaussian process learning, where the underlying Gaussian process is a discrete Gaussian Markov Random Field. This speed-up is accomplished by employing smart computational linear algebra, state-of-the-art algorithms, and a careful divide-and-conquer evaluation strategy. Problems of significantly greater size than any other existing algorithm with similar guarantees can solve are solved as illustrations.

Estimating portfolio risk measures and classifying portfolio risk levels in real time are important yet challenging tasks. In this paper, we propose to build a logistic regression model using data generated in past simulation experiments and to use the model to predict portfolio risk measures and classify risk levels at any time. We further explore regularization techniques, simulation model structure, and additional simulation budget to enhance the estimators of the logistic regression model to make its predictions more precise. Our numerical results show that the proposed methods work well. Our work may be viewed as an example of the recently proposed idea of simulation analytics, which treats a simulation model as a data generator and proposes to apply data analytics tools to the simulation outputs to uncover conditional statements. Our work shows that the simulation analytics idea is viable and promising in the field of financial risk management.

“Virtual statistics,” as we define them, are estimators of performance measures that are conditional on the occurrence of an event; virtual waiting time of a customer arriving to a queue at time [Formula: see text] is one example of virtual performance. In this paper, we describe a [Formula: see text]-nearest-neighbor method for estimating virtual performance postsimulation from the retained sample paths, examining both its small-sample and asymptotic properties and providing two approaches for measuring the error of the [Formula: see text]-nearest-neighbor estimator. We implement leave-one-replication-out cross-validation for tuning a single parameter [Formula: see text] to use for any time (or times) of interest and evaluate the prediction performance of the [Formula: see text]-nearest-neighbor estimator via controlled studies. As a by-product, this paper motivates a different way of thinking about how to process the output from dynamic, discrete-event simulation.

Selecting the solution with the largest or smallest mean of a primary performance measure from a finite set of solutions while requiring secondary performance measures to satisfy certain constraints is called constrained selection of the best (CSB) in the simulation ranking and selection literature. In this paper, we consider CSB problems with secondary performance measures that must satisfy probabilistic constraints, and we call such problems chance constrained selection of the best (CCSB). We design procedures that first check the feasibility of all solutions and then select the best among all the sample feasible solutions. We prove the statistical validity of these procedures for variations of the CCSB problem under the indifference-zone formulation. Numerical results show that the proposed procedures can efficiently handle CCSB problems with up to 100 solutions, each with five chance constraints.

We consider the problem of deriving confidence intervals for the mean response of a system that is represented by a stochastic simulation whose parametric input models have been estimated from “real-world” data. As opposed to standard simulation confidence intervals, we provide confidence intervals that account for uncertainty about the input model parameters; our method is appropriate when enough simulation effort can be expended to make simulation-estimation error relatively small. To achieve this we introduce metamodel-assisted bootstrapping that propagates input variability through to the simulation response via an equation-based model rather than by simulating. We develop a metamodel strategy and associated experiment design method that avoid the need for low-order approximation to the response and that minimizes the impact of intrinsic (simulation) error on confidence level accuracy. Asymptotic analysis and empirical tests over a wide range of simulation effort show that confidence intervals obtained via metamodel-assisted bootstrapping achieve the desired coverage.

We propose an adaptive hyperbox algorithm (AHA), which is an instance of a locally convergent, random search algorithm for solving discrete optimization via simulation problems. Compared to the COMPASS algorithm, AHA is more efficient in high-dimensional problems. By analyzing models of the behavior of COMPASS and AHA, we show why COMPASS slows down significantly as dimension increases, whereas AHA is less affected. Both AHA and COMPASS can be used as the local search algorithm within the Industrial Strength COMPASS framework, which consists of a global search phase, a local search phase, and a final cleanup phase. We compare the performance of AHA to COMPASS within the framework of Industrial Strength COMPASS and as stand-alone algorithms. Numerical experiments demonstrate that AHA scales up well in high-dimensional problems and has similar performance to COMPASS in low-dimensional problems.

Controlled sequential bifurcation (CSB) is a factor-screening method for discrete-event simulations. It combines a multistage hypothesis testing procedure with the original sequential bifurcation procedure to control both the power for detecting important effects at each bifurcation step and the Type I error for each unimportant factor under heterogeneous variance conditions when a main-effects model applies. This paper improves the CSB procedure in two aspects. First, a new fully sequential hypothesis-testing procedure is introduced that greatly improves the efficiency of CSB. Moreover, this paper proposes CSB-X, a more general CSB procedure that has the same error control for screening main effects that CSB does, even when two-factor interactions are present. The performance of the new method is proven and compared with the original CSB procedure.

Simulation models of real-life systems often assume stationary (homogeneous) Poisson arrivals. Therefore, when nonstationary arrival processes are required, it is natural to assume Poisson arrivals with a time-varying arrival rate. For many systems, however, this provides an inaccurate representation of the arrival process that is either more or less variable than Poisson. In this paper we extend techniques that transform a stationary Poisson arrival process into a nonstationary Poisson arrival process (NSPP) by transforming a stationary renewal process into a nonstationary, non-Poisson (NSNP) arrival process. We show that the desired arrival rate is achieved and that when the renewal base process is either more or less variable than Poisson, then the NSNP process is also more or less variable, respectively, than an NSPP. We also propose techniques for specifying the renewal base process when presented properties of, or data from, an arrival process and illustrate them by modeling real arrival data.

Time-series input processes occur naturally in the stochastic simulation of many service, communications, and manufacturing systems, and there are a variety of time-series input models available to match a given collection of properties, typically a marginal distribution and an autocorrelation structure specified via the use of one or more time lags. The focus of this paper is the situation in which the collection of properties are not “given,” but data are available from which a time-series input model is to be estimated. The input model we consider is the very flexible autoregressive-to-anything (ARTA) model of Cario and Nelson [Cario, M. C., B. L. Nelson. 1996. Autoregressive to anything: Time-series input processes for simulation. Oper. Res. Lett. 19 51–58]. Recently, we developed a statistically valid algorithm (ARTAFIT) for fitting this model to stationary univariate time-series data using marginal distributions from the Johnson translation system. In this paper, we perform a comprehensive numerical study to assess the performance of our algorithm relative to the two most commonly used approaches: (a) fitting the marginal distribution but ignoring the autocorrelation structure, and (b) fitting separately the marginal distribution as in (a) and the autocorrelation structure using the sample autocorrelation function. We find that ARTAFIT, which fits the marginal distribution and the autocorrelation structure jointly, outperforms both (a) and (b), and we demonstrate the importance of taking dependencies into account while developing input models for stochastic simulation.

Cycle time-throughput (CT-TH) percentile curves quantify the relationship between percentiles of cycle time and factory throughput, and they can play an important role in strategic planning for manufacturing systems. In this paper, a highly flexible distribution, the generalized gamma, is used to represent the underlying distribution of cycle time. To obtain CT-TH percentile curves, we use a factory simulation to fit metamodels for the first three CT-TH moment curves throughout the throughput range of interest, determine the parameters of the generalized gamma by matching moments, and obtain any percentile of interest by inverting the distribution. To insure efficiency and control estimation error, simulation experiments are built up sequentially using a multistage procedure. Numerical results are presented to demonstrate the effectiveness of the approach.

We develop a numerically exact method for evaluating the time-dependent mean, variance, and higher-order moments of the number of entities in a Pht/Pht/∞ queueing system. We also develop a numerically exact method for evaluating the distribution function and moments of the virtual sojourn time for any time t; in our setting, the virtual sojourn time is equivalent to the service time for virtual entities arriving to the system at that time t. We include several examples using software that we have developed and have put in downloadable form in the Online Supplement to this paper on the journal's website.

We demonstrate a numerically exact method for evaluating the time-dependent mean, variance, and higher-order moments of the number of entities in the multiclass [Ph t /Ph t /∞] K queueing network system, as well as at the individual network nodes. We allow for multiple, independent, time-dependent entity classes and develop time-dependent performance measures by entity class at the nodal and network levels. We also demonstrate a numerically exact method for evaluating the distribution function and moments of virtual sojourn time through the network for virtual entities, by entity class, arriving to the system at time t. We include an example using software that we have developed and have put in downloadable form in the Online Supplement to this paper on the journal's website.

We present and evaluate three ranking-and-selection procedures for use in steady-state simulation experiments when the goal is to find which among a finite number of alternative systems has the largest or smallest long-run average performance. All three procedures extend existing methods for independent and identically normally distributed observations to general stationary output processes, and all procedures are sequential. We also provide our thoughts about the evaluation of simulation design and analysis procedures, and illustrate these concepts in our evaluation of the new procedures.

An ARTA (AutoRegressive-to-Anything) Process is a time series with arbitrary marginal distribution and autocorrelation structure specified through finite lag p. We develop an efficient numerical method for fitting ARTA processes and discuss its implementation in the software ARTAFACTS. We also present the software ARTAGEN that generates observations from ARTA processes for use as inputs to a computer simulation. We illustrate the use of the software with a real-world example.