WSC'00
Final Program Abstracts
Analysis Methodology I Track

Using Lattice Rules for Variance Reduction in Simulation
Christiane Lemieux (University of Calgary) and Pierre L'Ecuyer (Université de Montréal)

Abstract:
Quasi-Monte Carlo methods are designed to improve upon the Monte Carlo method for multidimensional numerical integration by using a more regularly distributed point set than the i.i.d. sample associated with Monte Carlo. Lattice rules are one family of quasi-Monte Carlo methods, originally proposed by Korobov in 1959. In this paper, we explain how randomized lattice rules can be used to construct efficient estimators for typical simulation problems, and we give several numerical examples. We are interested in two main aspects: Studying the variance of these estimators and finding which properties of the lattice rules should be considered when defining a selection criterion to rate and choose them. Our numerical results for three different problems illustrate how this methodology typically improves upon the usual Monte Carlo simulation method.

Quasi-Monte Carlo Methods in Cash Flow Testing Simulations
Michael G. Hilgers (University of Missouri - Rolla)

Abstract:
What actuaries call cash flow testing is a large-scale simulation pitting a company's current policy obligation against future earnings based on interest rates. While life contingency issues associated with contract payoff are a mainstay of the actuarial sciences, modeling the random fluctuations of US Treasury rates is less studied. Furthermore, applying standard simulation techniques, such as the Monte Carlo method, to actual multi-billion dollar companies produce a simulation that can be computationally prohibitive. In practice, only hundreds of sample paths can be considered, not the usual hundreds of thousands one might expect for a simulation of this complexity. Hence, insurance companies have a desire to accelerate the convergence of the estimation procedure. This paper reports the results of cash flow testing simulations performed for Conseco L.L.C. using so-called quasi-Monte Carlo techniques. In these, pseudo-random number generation is replaced with deterministic low discrepancy sequences. It was found that by judicious choice of subsequences, the quasi-Monte Carlo method provided a consistently tighter estimate, than the traditional methods, for a fixed, small number of sample paths. The techniques used to select these subsequences are discussed.

Generating 'Dependent' Quasi-Random Numbers
Shane G. Henderson (University of Michigan) and Belinda A. Chiera and Roger M. Cooke (Delft University of Technology)

Abstract:
Under certain conditions on the integrand, quasi-Monte Carlo methods for estimating integrals (expectations) converge faster asymptotically than Monte Carlo methods. Motivated by this result we consider the generation of quasi-random vectors with given marginals and given correlation matrix. We extend the "Normal To Anything" (NORTA) method, introduced by Cario and Nelson, to this context, and term the extension the "Quasi-Random to Anything" (QUARTA) method.

Adaptively Choosing the Best Procedure for Selecting the Best System
Justin Boesel (The MITRE Corporation)

Abstract:
This paper presents a method that uses initial sample data to choose between statistical procedures for identifying the simulated system with the best (maximum or minimum) expected performance. The method chooses the procedure that minimizes the additional number of simulation replications required to return a pre-specified probability guarantee. This problem may be encountered after a heuristic search procedure has been applied in a simulation-optimization context. In this setting, initial samples from each system may already have been taken, but because of stochastic variation, the system with the best sample mean at the end of the search procedure may not be the true best system encountered during the search. Empirical work in previous papers suggests that the relative number of additional replications required by existing procedures depends on factors --- such as the configuration of the systems' means and their variances --- that may be unknown prior to initial data collection. These results motivated the approach taken in this paper, where we postpone the choice between statistical procedures until after observing the initial data.

Ranking and Selection for Steady-State Simulation
David Goldsman (Georgia Institute of Technology), William S. Marshall (Georgia Tech Research Institute) and Seong-Hee Kim and Barry L. Nelson (Northwestern University)

Abstract:
We present and evaluate two 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. Both procedures extend existing methods for independent and identically normally distributed observations to general stationary output processes, and both procedures are sequential.

New Results on Procedures that Select the Best System Using CRN
Stephen E. Chick and Koichiro Inoue (The University of Michigan)

Abstract:
This tutorial describes some ways that Bayesian methods address problems that arise during simulation studies. This includes quantifying uncertainty about input distributions and parameters, sensitivity analysis, and the selection of the best of several simulated alternatives. Focus is on illustrating the main ideas and their relevance to practical problems. Numerous citations for both introductory and more advanced material provide a launching pad into the Bayesian literature.

Nonparametric Estimation of Nonhomogeneous Poisson Processes Using Wavelets
Michael E. Kuhl and Prashant S. Bhairgond (Louisaina State University)

Abstract:
Nonhomogeneous Poisson processes (NHPPs) are frequently used in stochastic simulations to model nonstationary point processes. These NHPP models are often constructed by estimating the rate function from one or more observed realizations of the process. Both parametric and nonparametric models have been developed for the NHPP rate function. The current parametric models require prior knowledge of the behavior of the NHPP under study for model selection. The current nonparametric estimators, in general, require the storage of all of the observed data. Other hybrid approaches have also been developed. This paper focuses on the nonparametric estimation of the rate function of a nonhomogeneous Poisson process using wavelets. The advantages of wavelets include the flexibil-ity of a nonparametric estimator enabling one to model the nonstationary rate function of an NHPP without prior knowledge or assumptions about the behavior of the proc-ess. Furthermore, this method has some advantages of current nonparametric techniques. Thus, using wavelets we can develop an efficient yet highly flexible NHPP rate function. In this paper, we develop the methodology required for constructing a wavelet estimator for the NHPP rate function. In addition, we present an experimental performance evaluation for this method. Nonhomogeneous Poisson processes (NHPPs) are frequently used in stochastic simulations to model nonstationary point processes. These NHPP models are often constructed by estimating the rate function from one or more observed realizations of the process. Both parametric and nonparametric models have been developed for the NHPP rate function. The current parametric models require prior knowledge of the behavior of the NHPP under study for model selection. The current nonparametric estimators, in general, require the storage of all of the observed data. Other hybrid approaches have also been developed. This paper focuses on the nonparametric estimation of the rate function of a nonhomogeneous Poisson process using wavelets. The advantages of wavelets include the flexibility of a nonparametric estimator enabling one to model the nonstationary rate function of an NHPP without prior knowledge or assumptions about the behavior of the process. Furthermore, this method has some advantages of current nonparametric techniques. Thus, using wavelets we can develop an efficient yet highly flexible NHPP rate function. In this paper, we develop the methodology required for constructing a wavelet estimator for the NHPP rate function. In addition, we present an experimental performance evaluation for this method.

Simulation from Non-Standard Distributions Using Envelope Methods
Michael J. Evans (University of Toronto) and Tim B. Swartz (Simon Fraser University)

Abstract:
This paper considers the development of envelope methods as a tool for simulation. Envelope methods are based on the construction of simple envelopes to functions. The proposed envelopes are general, require little input from the user and are based on the concavity structure of the function or some transformation of the function. The construction of these envelopes facilitates variate generation using the adaptive rejection algorithm.

Input Modeling Using a Computer Algebra System
Diane L. Evans and Lawrence M. Leemis (College of William and Mary)

Abstract:
Input modeling that involves fitting standard univariate parametric probability distributions is typically performed using an input modeling package. These packages typically fit several distributions to a data set, then determine the distribution with the best fit by comparing goodness-of-fit statistics. But what if an appropriate input model is not included in one of these packages? The modeler must resort to deriving the appropriate estimators by hand for the appropriate input model. The purpose of this paper is to investigate the use of a prototype Maple-based probability language, known as APPL (A Probability Programming Language), for input modeling. This language allows an analyst to specify a standard or non-standard distribution for an input model, and have the derivations performed automatically. Input modeling serves as an excellent arena for illustrating the applicability and usefulness of APPL. Besides including pre-defined types for over 45 different continuous and discrete random variables and over 30 procedures for manipulating random variables (e.g., convolution, transformation), APPL contains input modeling procedures for parameter estimation, plotting empirical and fitted CDFs, and performing goodness-of-fit tests. Using examples, we illustrate its utility for input modeling.

Importance Sampling in Derivative Securities Pricing
Yi Su and Michael C. Fu (University of Maryland)

Abstract:
We formulate the importance sampling problem as a parametric minimization problem under the original measure and use a combination of infinitesimal perturbation analysis (IPA) and stochastic approximation (SA) to minimize the variance of the price estimation. Compared to existing methods, the IPA estimator derived in this paper has significantly smaller estimation variance and doesn't depend on the form of payoff functions and differentiability of the sample path, and thus is more universally applicable and computationally efficient. Under suitable conditions, the objective function is a convex function, the IPA estimator presented is unbiased, and the corresponding stochastic approximation algorithm converges to the true optimal value.

A Real Options Design for Quality Control Charts
Harriet Black Nembhard, Leyuan Shi, and Mehmet Aktan (University of Wisconsin-Madison)

Abstract:
We develop a financial model for a manufacturing process where quality can be affected by an assignable cause. We evaluate the options associated with applying a statistical process control chart using pentanomial lattice and Monte Carlo simulation methods. By connecting the aspects of market dynamics with the manufacturing operational aspects, we now have a way to help decision makers address the bottom-line profitability associated with the quality control decision.

Variance Reduction Techniques for Value-at-Risk with Heavy-Tailed Risk Factors
Paul Glasserman (Columbia University), Philip Heidelberger (IBM Research Division) and Perwez Shahabuddin (Columbia University)

Abstract:
The calculation of value-at-risk (VAR) for large portfolios of complex instruments is among the most demanding and widespread computational challenges facing the financial industry. Current methods for calculating VAR include comparatively fast numerical approximations---especially the linear and quadratic (delta-gamma) approximations---and more robust but more computationally demanding Monte Carlo simulation. The linear and delta-gamma methods typically rely on an assumption that the underlying market risk factors have a Gaussian distribution over the VAR horizon. But there is ample empirical evidence that market data is more accurately described by heavy-tailed distributions. Capturing heavy tails in VAR calculations has to date required highly time-consuming Monte Carlo simulation. We describe two methods for computationally efficient calculation of VAR in the presence of heavy-tailed risk factors, specifically when risk factors have a multivariate t distribution. The first method uses transform inversion to develop a fast numerical algorithm for computing the distribution of the heavy-tailed delta-gamma approximation. For greater accuracy, the second method uses the numerical approximation to guide in the design of an effective Monte Carlo variance reduction technique; the algorithm combines importance sampling and stratified sampling. This method can produce enormous speed-ups compared with standard Monte Carlo.

Integrating Optimization and Simulation: Research and Practice
Michael C. Fu (University of Maryland), Sigrún Andradóttir (Georgia Institute of Technology), John S. Carson (AutoSimulations, Inc.), Fred Glover (University of Colorado), Charles R. Harrell (Brigham Young University), Yu-Chi Ho (Harvard University), James P. Kelly (University of Colorado) and Stephen M. Robinson (University of Wisconsin at Madison)

Abstract:
We formulate the importance sampling problem as a parametric minimization problem under the original measure and use a combination of infinitesimal perturbation analysis (IPA) and stochastic approximation (SA) to minimize the variance of the price estimation. Compared to existing methods, the IPA estimator derived in this paper has significantly smaller estimation variance and doesn't depend on the form of payoff functions and differentiability of the sample path, and thus is more universally applicable and computationally efficient. Under suitable conditions, the objective function is a convex function, the IPA estimator presented is unbiased, and the corresponding stochastic approximation algorithm converges to the true optimal value.

A Stopping Procedure Based on Phi-Mixing Conditions
E. Jack Chen and W. David Kelton (University of Cincinnati)

Abstract:
The use of batch means is a well-known technique for estimating the variance of mean point estimators computed from a simulation experiment. This paper discusses implementation of a sequential procedure to determine the batch size for constructing confidence intervals for a simulation estimator of the steady-state mean of a stochastic process. Our quasi-independent (QI) procedure increases the batch size progessively until a certain number of essentially i.i.d. samples are obtained. We show that our sequential procedure gives valid confidence intervals. The only (mild) assumption is that the correlations of the stochastic process output sequence die off. An experimental performance evaluation demonstrates the validity of the QI procedure.

Experimental Performance Evaluation of Batch Means Procedures for Simulation Output Analysis
Natalie M. Steiger (University of Maine) and James R. Wilson (North Carolina State University)

Abstract:
We summarize the results of an extensive experimental performance evaluation of selected batch means procedures for building a confidence interval for a steady-state expected simulation response. We compare the performance of the well-known ABATCH and LBATCH procedures versus ASAP, a recently proposed variant of the method of nonoverlapping batch means (NOBM) that operates as follows: the batch size is progressively increased until either (a) the batch means pass the von Neumann test for independence, and then ASAP delivers a classical NOBM confidence interval; or (b) the batch means pass the Shapiro-Wilk test for multivariate normality, and then ASAP delivers a correlation-adjusted confidence interval. The latter correction is based on an inverted Cornish-Fisher expansion for the classical NOBM t-ratio, where the terms of the expansion are estimated via an autoregressive--moving average time series model of the batch means. Applying ABATCH, ASAP, and LBATCH to the analysis of a suite of twenty test problems involving discrete-time Markov chains, time-series processes, and queueing systems, we found ASAP to deliver confidence intervals that not only satisfy a user-specified absolute or relative precision requirement but also frequently outperform the corresponding confidence intervals delivered by ABATCH and LBATCH with respect to coverage probability.

Simulation Output Analysis via Dynamic Batch Means
Yingchieh Yeh and Bruce Schmeiser (Purdue University)

Abstract:
This paper is focused on estimating the quality of the sample mean from a steady-state simulation experiment with consideration of computational efficiency, memory requirement, and statistical efficiency. In addition, we seek methods that do not require knowing run length a priori. We develop an algorithm of nonoverlapping batch means that is implemented in fixed memory by dynamically changing both batch size and number of batches as the simulation runs. The algorithm, denoted by DBM for Dynamic Batch Means, requires computation time similar to other batch means data-collection methods, despite its fixed memory requirement. To achieve satisfactory statistical efficiency of DBM, we propose two associated estimators, V_TBM and V_PBM, of the variance of the sample mean and investigate their statistical properties. Our study shows that the estimator V_PBM with parameter w=1 is, as a practical matter, better than the other.

Adaptive Importance Sampling Simulation of Queueing Networks
Pieter-Tjerk de Boer and Victor F. Nicola (University of Twente) and Reuven Y. Rubinstein (Technion)

Abstract:
In this paper, a method is presented for the efficient estimation of rare-event (overflow) probabilities in Jackson queueing networks using importance sampling. The method differs in two ways from methods discussed in most earlier literature: the change of measure is state-dependent, i.e., it is a function of the content of the buffers, and the change of measure is determined using a cross-entropy-based adaptive procedure. This method yields asymptotically efficient estimation of overflow probabilities of queueing models for which it has been shown that methods using a state-independent change of measure are not asymptotically efficient. Numerical results demonstrating the effectiveness of the method are presented as well.

Simulating GI/GI/1 Queues and Insurance Risk Processes with Subexponential Distributions
Nam Kyoo Boots (Vrije Universiteit) and Perwez Shahabuddin (Columbia University)

Abstract:
This paper deals with estimating small tail probabilities of the steady-state waiting time in a GI/GI/1 queue with heavy-tailed (subexponential) service times. The problem of estimating infinite horizon ruin probabilities in insurance risk processes with heavy-tailed claims can be transformed into the same framework. It is well-known that naive simulation is ineffective for estimating small probabilities and special fast simulation techniques like importance sampling, multilevel splitting, etc., have to be used. Previous fast simulation techniques for queues with subexponential service times have been confined to M/GI/1 queueing systems. The general approach is to use the Pollaczek-Khintchine transformation to transform the problem into that of estimating the tail distribution of a geometric sum of independent subexponential random variables. However, no such useful transformation exists when one goes from Poisson arrivals to general interarrival-time distributions. We describe an approach that is based on directly simulating the random walk associated with the waiting-time process of the GI/GI/1 queue, using a change of measure called delayed subexponential twisting -- an importance sampling idea recently developed and found useful in the context of M/GI/1 heavy-tailed simulations.

Practical Aspects of Simulating Systems having Arrival Processes with Long-Range Dependence
Robert Geist and James Westall (Clemson University)

Abstract:
Analysis of network traffic indicates that packet arrival processes have significant stochastic dependence. It has been suggested that this dependence is so strong as to be well-modeled by long-range dependent processes. The fact that no finite process can be said to be truly long-range dependent poses potentially serious obstacles to simulation modeling. In this paper, we explore some of these obstacles and propose practical methods for obtaining useful results with simulations of manageable duration. Keywords: Analysis methodology; network traffic models; stochastic dependence; long-range dependence.

Automatic Random Variate Generation for Simulation Input
W. Hörmann (Bogaziçi University Istanbul) and J. Leydold (University of Economics and BA Vienna)

Abstract:
We develop and evaluate algorithms for generating random variates for simulation input. One group called automatic, or black-box algorithms can be used to sample from distributions with known density. They are based on the rejection principle. The hat function is generated automatically in a setup step using the idea of transformed density rejection. There the density is transformed into a concave function and the minimum of several tangents is used to construct the hat function. The resulting algorithms are not too complicated and are quite fast. The principle is also applicable to random vectors. A second group of algorithms is presented that generate random variates directly from a given sample by implicitly estimating the unknown distribution. The best of these algorithms are based on the idea of naive resampling plus added noise. These algorithms can be interpreted as sampling from the kernel density estimates. This method can be also applied to random vectors. There it can be interpreted as a mixture of naive resampling and sampling from the multi-normal distribution that has the same covariance matrix as the data. The algorithms described in this paper have been implemented in ANSI C in a library called UNURAN which is available via anonymous ftp.

Fast Combined Multiple Recursive Generators with Multipliers of the Form a = ± 2^q ± 2^r
Pierre L'Ecuyer and Renée Touzin (Université de Montréal)

Abstract:
We study a class of combined multiple recursive random number generators constructed in a way that each component runs fast and is easy to implement, while the combination enjoys excellent structural properties as measured by the spectral test. Each component is a linear recurrence of order k > 1, modulo a large prime number, and the coefficients are either 0 or are of the form a = ± 2^q or a = ± 2^q ± 2^r. This allows a simple and very fast implementation, because each modular multiplication by a power of 2 can be implemented via a shift, plus a few additional operations for the modular reduction. We select the parameters in terms of the performance of the combined generator in the spectral test. We provide a specific implementation.

A New Class of Linear Feedback Shift Register Generators
Pierre L'Ecuyer and Francois Panneton (Université de Montréal)

Abstract:
An efficient implementation of linear feedback shift register sequences with a given characteristic polynomial is obtained by a new method. It involves a polynomial linear congruential generator over the finite field with two elements. We obtain maximal equidistribution by constructing a suitable output mapping. Local randomness could be improved by combining the generator's output with that of some other (e.g., nonlinear and efficient) generator.

Analysis of Simulation Factorial Experiments by EDF Resample Statistics
Russell C. H. Cheng and Owen D. Jones (University of Southampton)

Abstract:
The output from simulation factorial experiments can be complex and may not be amenable to standard methods of estimation like ANOVA. We consider the situation where the simulation output may not satisfy normality assumptions, but more importantly, where there may be differences in output at different factor combinations, but these are not simply differences in means. We show that EDF statistics can provide a similar but potentially more sensitive analysis to that provided by ANOVA. Moreover we show that with the use of resampling, we can generate accurate critical values for tests of hypothesis under much weaker conditions than those required for ANOVA tests. The method is illustrated with an example based on an actual simulation experiment comparing two methods of operating a production facility under different production levels.

Low Cost Response Surface Methods for and from Simulation Optimization
Theodore Allen (The Ohio State University) and Liyang Yu (i2 Technologies, Inc.)

Abstract:
We propose "low cost response surface methods" (LCRSM) that typically require half the experimental runs of standard response surface methods based on central composite and Box Behnken designs but yield comparable or lower modeling errors under realistic assumptions. In addition, the LCRSM methods have substantially lower modeling errors and greater expected savings compared with alternatives with comparable numbers of runs, including small composite designs and computer-generated designs based on popular criteria such as D-optimality. Therefore, when simulation runs are expensive, low cost response surface methods can be used to create regression meta-models for queuing or other system optimization. The LCRSM procedures are also apparently the first experimental design methods derived as the solution to a simulation optimization problem. For these reasons, we say that LCRSM are "for and from" simulation optimization. We compare the proposed LCRSM methods with a large number of alternatives based on six criteria. We conclude that the proposed methods offer attractive alternative to current methods in many relevant situations.

Simulation Optimization of Stochastic Systems with Integer Variables by Sequential Linearization
S. J. Abspoel, L. F. P. Etman, J. Vervoort, and J. E. Rooda (Eindhoven University of Technology)

Abstract:
Discrete-event simulation is widely used to analyse and improve the performance of manufacturing systems. The related optimization problem often includes integer design variables and is defined by objective function and constraints that are expected values of stochastic functions. These stochastic functions have to be evaluated via the simulation model at the discrete levels of the integer design parameters. For such a simulation optimization problem with integer variables we have developed an optimization strategy that is based on a series of linear approximate subproblems. Each subproblem is built from the outcomes of simulation experiments. A D-optimal design of experiments is used to plan the simulation experiments. Stochasticity in constraint and objective functions is dealt with explicitly using safety indices. Two test problems will be presented to illustrate the optimization strategy. This includes a simulation based four-station production flow line problem.