We present a numerical method to compute the survival function and the moments of the exit time for a piecewise-deterministic Markov process (PDMP). Our approach is based on the quantization of an underlying discrete-time Markov chain related to the PDMP. The approximation we propose is easily computable and is even flexible with respect to the exit time we consider. We prove the convergence of the algorithm and obtain bounds for the rate of convergence in the case of the moments. We give an academic example and a model from the reliability field to illustrate the results of the paper.

We derive an expansion for the (expected) difference between the continuously monitored supremum and evenly monitored discrete maximum over a finite time horizon of a jump diffusion process with independent and identically distributed normal jump sizes. The monitoring error is of the form a 0 / N 1/2 + a 1 / N 3/2 + · · · + b 1 / N + b 2 / N 2 + b 4 / N 4 + · · ·, where N is the number of monitoring intervals. We obtain explicit expressions for the coefficients {a 0, a 1, …, b 1, b 2, …}. In particular, a 0 is proportional to the value of the Riemann zeta function at ½, a well-known fact that has been observed for Brownian motion in applied probability and mathematical finance.

We propose a method for estimating first passage time densities of one-dimensional diffusions via Monte Carlo simulation. Our approach involves a representation of the first passage time density as the expectation of a functional of the three-dimensional Brownian bridge. As the latter process can be simulated exactly, our method leads to almost unbiased estimators. Furthermore, since the density is estimated directly, a convergence of order 1 / √N, where N is the sample size, is achieved, which is in sharp contrast to the slower nonparametric rates achieved by kernel smoothing of cumulative distribution functions.

The variance minimization (VM) and cross-entropy (CE) methods are two versatile adaptive importance sampling procedures that have been successfully applied to a wide variety of difficult rare-event estimation problems. We compare these two methods via various examples where the optimal VM and CE importance densities can be obtained analytically. We find that in the cases studied both VM and CE methods prescribe the same importance sampling parameters, suggesting that the criterion of minimizing the CE distance is very close, if not asymptotically identical, to minimizing the variance of the associated importance sampling estimator.

We present an exact simulation algorithm for the stationary distribution of the customer delay D for first-in–first-out (FIFO) M/G/c queues in which ρ=λ/μ<1. We assume that the service time distribution G(x)=P(S≤x),x≥0 (with mean 0<E(S)=1/μ<∞), and its corresponding equilibrium distribution Ge(x)=μ∫0x P(S>y)dy are such that samples of them can be simulated. We further assume that G has a finite second moment. Our method involves the general method of dominated coupling from the past (DCFTP) and we use the single-server M/G/1 queue operating under the processor sharing discipline as an upper bound. Our algorithm yields the stationary distribution of the entire Kiefer–Wolfowitz workload process, the first coordinate of which is D. Extensions of the method to handle simulating generalized Jackson networks in stationarity are also remarked upon.

A stochastic perpetuity takes the form D∞=∑n=0∞ exp(Y1+⋯+Yn)Bn, where Yn:n≥0) and (Bn:n≥0) are two independent sequences of independent and identically distributed random variables (RVs). This is an expression for the stationary distribution of the Markov chain defined recursively by Dn+1=AnDn+Bn, n≥0, where An=eYn; D∞ then satisfies the stochastic fixed-point equation D∞D̳AD∞+B, where A and B are independent copies of the An and Bn (and independent of D∞ on the right-hand side). In our framework, the quantity Bn, which represents a random reward at time n, is assumed to be positive, unbounded with EBnp <∞ for some p>0, and have a suitably regular continuous positive density. The quantity Yn is assumed to be light tailed and represents a discount rate from time n to n-1. The RV D∞ then represents the net present value, in a stochastic economic environment, of an infinite stream of stochastic rewards. We provide an exact simulation algorithm for generating samples of D∞. Our method is a variation of dominated coupling from the past and it involves constructing a sequence of dominating processes.

We consider the problem of estimating the probability that the maximum of a Gaussian process with negative mean and indexed by positive integers reaches a high level, say b. In great generality such a probability converges to 0 exponentially fast in a power of b. Under mild assumptions on the marginal distributions of the process and no assumption on the correlation structure, we develop an importance sampling procedure, called the target bridge sampler (TBS), which takes a polynomial (in b) number of function evaluations to achieve a small relative error. The procedure also yields samples of the underlying process conditioned on hitting b in finite time. In addition, we apply our method to the problem of estimating the tail of the maximum of a superposition of a large number, n, of independent Gaussian sources. In this situation TBS achieves a prescribed relative error with a bounded number of function evaluations as n ↗ ∞. A remarkable feature of TBS is that it is not based on exponential changes of measure. Our numerical experiments validate the performance indicated by our theoretical findings.

In a fork-join network each incoming job is split into K tasks and the K tasks are simultaneously assigned to K parallel service stations for processing. For the distributions of response times and queue lengths of fork-join networks, no explicit formulae are available. Existing methods provide only analytic approximations for the response time and the queue length distributions. The accuracy of such approximations may be difficult to justify for some complicated fork-join networks. In this paper we propose a perfect simulation method based on coupling from the past to generate exact realisations from the equilibrium of fork-join networks. Using the simulated realisations, Monte Carlo estimates for the distributions of response times and queue lengths of fork-join networks are obtained. Comparisons of Monte Carlo estimates and theoretical approximations are also provided. The efficiency of the sampling algorithm is shown theoretically and via simulation.

Mixed-level orthogonal arrays are basic structures in experimental design. We develop three algorithms that compute Rao- and Gilbert-Varshamov-type bounds for mixed-level orthogonal arrays. The computational complexity of the terms involved in the original combinatorial representations of these bounds can grow fast as the parameters of the arrays increase and this justifies the construction of these algorithms. The first is a recursive algorithm that computes the bounds exactly, the second is based on an asymptotic analysis, and the third is a simulation algorithm. They are all based on the representation of the combinatorial expressions that appear in the bounds as expectations involving a symmetric random walk. The Markov property of the underlying random walk gives the recursive formula to compute the expectations. A large deviation (LD) analysis of the expectations provides the asymptotic algorithm. The asymptotically optimal importance sampling (IS) of the same expectation provides the simulation algorithm. Both the LD analysis and the construction of the IS algorithm use a representation of these problems as a sequence of stochastic optimal control problems converging to a limit calculus of a variations problem. The construction of the IS algorithm uses a recently discovered method of using subsolutions to the Hamilton-Jacobi-Bellman equations associated with the limit problem.

Max-stable random fields play a central role in modeling extreme value phenomena. We obtain an explicit formula for the conditional probability in general max-linear models, which include a large class of max-stable random fields. As a consequence, we develop an algorithm for efficient and exact sampling from the conditional distributions. Our method provides a computational solution to the prediction problem for spectrally discrete max-stable random fields. This work offers new tools and a new perspective to many statistical inference problems for spatial extremes, arising, for example, in meteorology, geology, and environmental applications.

We consider a feed-forward network with a single-server station serving jobs with multiple levels of priority. The service discipline is preemptive in that the server always serves a job with the current highest level of priority. For this system with discontinuous dynamics, we establish the sample path large deviation principle using a weak convergence argument. In the special case where jobs have two different levels of priority, we also explicitly identify the exponential decay rate of the total population overflow probabilities by examining the geometry of the zero-level sets of the system Hamiltonians.

We present a method for computing the probability density function (PDF) and the cumulative distribution function (CDF) of a nonnegative infinitely divisible random variable X. Our method uses the Lévy-Khintchine representation of the Laplace transform Ee-λX = e-ϕ(λ), where ϕ is the Laplace exponent. We apply the Post-Widder method for Laplace transform inversion combined with a sequence convergence accelerator to obtain accurate results. We demonstrate this technique on several examples, including the stable distribution, mixtures thereof, and integrals with respect to nonnegative Lévy processes.

In this paper we prove that the stationary distribution of populations in genetic algorithms focuses on the uniform population with the highest fitness value as the selective pressure goes to ∞ and the mutation probability goes to 0. The obtained sufficient condition is based on the work of Albuquerque and Mazza (2000), who, following Cerf (1998), applied the large deviation principle approach (Freidlin-Wentzell theory) to the Markov chain of genetic algorithms. The sufficient condition is more general than that of Albuquerque and Mazza, and covers a set of parameters which were not found by Cerf.

This paper demonstrates the application of a new higher-order weak approximation, called the Kusuoka approximation, with discrete random variables to non-commutative multi-factor models. Our experiments show that using the Heath–Jarrow–Morton model to price interest-rate derivatives can be practically feasible if the Kusuoka approximation is used along with the tree-based branching algorithm.

Geometric convergence to 0 of the probability that the goal has not been encountered by the nth generation is established for a class of genetic algorithms. These algorithms employ a quickly decreasing mutation rate and a crossover which restarts the algorithm in a controlled way depending on the current population and restricts execution of this crossover to occasions when progress of the algorithm is too slow. It is shown that without the crossover studied here, which amounts to a tempered restart of the algorithm, the asserted geometric convergence need not hold.

In a repulsive point process, points act as if they are repelling one another, leading to underdispersed configurations when compared to a standard Poisson point process. Such models are useful when competition for resources exists, as in the locations of towns and trees. Bertil Matérn introduced three models for repulsive point processes, referred to as types I, II, and III. Matérn used types I and II, and regarded type III as intractable. In this paper an algorithm is developed that allows for arbitrarily accurate approximation of the likelihood for data modeled by the Matérn type-III process. This method relies on a perfect simulation method that is shown to be fast in practice, generating samples in time that grows nearly linearly in the intensity parameter of the model, while the running times for more naive methods grow exponentially.

The waste-recycling Monte Carlo (WRMC) algorithm introduced by physicists is a modification of the (multi-proposal) Metropolis–Hastings algorithm, which makes use of all the proposals in the empirical mean, whereas the standard (multi-proposal) Metropolis–Hastings algorithm uses only the accepted proposals. In this paper we extend the WRMC algorithm to a general control variate technique and exhibit the optimal choice of the control variate in terms of the asymptotic variance. We also give an example which shows that, in contradiction to the intuition of physicists, the WRMC algorithm can have an asymptotic variance larger than that of the Metropolis–Hastings algorithm. However, in the particular case of the Metropolis–Hastings algorithm called the Boltzmann algorithm, we prove that the WRMC algorithm is asymptotically better than the Metropolis–Hastings algorithm. This last property is also true for the multi-proposal Metropolis–Hastings algorithm. In this last framework we consider a linear parametric generalization of WRMC, and we propose an estimator of the explicit optimal parameter using the proposals.

In this paper we study efficient simulation algorithms for estimating P(X›x), where X is the total time of a job with ideal time T that needs to be restarted after a failure. The main tool is importance sampling, where a good importance distribution is identified via an asymptotic description of the conditional distribution of T given X›x. If T≡t is constant, the problem reduces to the efficient simulation of geometric sums, and a standard algorithm involving a Cramér-type root, γ(t), is available. However, we also discuss an algorithm that avoids finding the root. If T is random, particular attention is given to T having either a gamma-like tail or a regularly varying tail, and to failures at Poisson times. Different types of conditional limit occur, in particular exponentially tilted Gumbel distributions and Pareto distributions. The algorithms based upon importance distributions for T using these asymptotic descriptions have bounded relative error as x→∞ when combined with the ideas used for a fixed t. Nevertheless, we give examples of algorithms carefully designed to enjoy bounded relative error that may provide little or no asymptotic improvement over crude Monte Carlo simulation when the computational effort is taken into account. To resolve this problem, an alternative algorithm using two-sided Lundberg bounds is suggested.

The standard Markov chain Monte Carlo method of estimating an expected value is to generate a Markov chain which converges to the target distribution and then compute correlated sample averages. In many applications the quantity of interest θ is represented as a product of expected values, θ = µ1 ⋯ µk, and a natural estimator is a product of averages. To increase the confidence level, we can compute a median of independent runs. The goal of this paper is to analyze such an estimator , i.e. an estimator which is a ‘median of products of averages’ (MPA). Sufficient conditions are given for to have fixed relative precision at a given level of confidence, that is, to satisfy . Our main tool is a new bound on the mean-square error, valid also for nonreversible Markov chains on a finite state space.