Quasi-Random Numbers and their Applications

Chair: David Humphrey (Nortel Networks)

*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.

Ranking, Selection, and Multiple Comparisons in Simulatioin

Chair: Daniel Ockerman (Manhattan Associates)

*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.

New Frontiers in Input Modeling

Chair: Paul Sanchez (Naval Postgraduate School)

*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.

Simulation in Financial Engineering

Chair: Laurel Travis (Metropolitan State University)

*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.

Panel on Integrating Optimization and Simulation: Research and Practice

Chair: Michael Fu (University of Maryland)

*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.

Batching Methods for Simulation Output Analysis

Chair: Bruce Schmeiser (Purdue University)

*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

Techniques for Simulating Difficult Queueing Problems

Chair: Michael R. Taaffe (University of Minnesota)

*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.

Random-Number and Random-Variate Generation

Chair: Saralees Nadarajah (University of California Santa Barbara)

*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.

Methods for Comparing and Optimizing Simulated Systems

Chair: Justin Boesel (The MITRE Corporation)

*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.