We introduce exact methods for the simulation of sample paths of one-dimensional diffusions with a discontinuity in the drift function. Our procedures require the simulation of finite-dimensional candidate draws from probability laws related to those of Brownian motion and its local time, and are based on the principle of retrospective rejection sampling. A simple illustration is provided.

Upon a set of backward orthogonal polynomials, we propose a novel multi-step numerical scheme for solving the decoupled forward-backward stochastic differential equations (FBSDEs). Under Lipschtiz conditions on the coefficients of the FBSDEs, we first get a general error estimate result which implies zero-stability of the proposed scheme, and then we further prove that the convergence rate of the scheme can be of high order for Markovian FBSDEs. Some numerical experiments are presented to demonstrate the accuracy of the proposed multi-step scheme and to numerically verify the theoretical results.

This paper concerns the computation of nonlinear crest distributions for irregular Stokes waves, and a numerical algorithm based on the Fast Fourier Transform (FFT) technique has been developed for carrying out the nonlinear computations. In order to further improve the computational efficiency, a new Transformed Rayleigh procedure is first proposed as another alternative for computing the nonlinear wave crest height distributions, and the corresponding computer code has also been developed. In the proposed Transformed Rayleigh procedure, the transformation model is chosen to be a monotonic exponential function, calibrated such that the first three moments of the transformed model match the moments of the true process. The numerical algorithm based on the FFT technique and the proposed Transformed Rayleigh procedure have been applied for calculating the wave crest distributions of a sea state with a Bretschneider spectrum and a sea statewith the surface elevation datameasured at the Poseidon platform. It is demonstrated in these two cases that the numerical algorithm based on the FFT technique and the proposed Transformed Rayleigh procedure can offer better predictions than those from using the empirical wave crest distribution models. Meanwhile, it is found that our proposed Transformed Rayleigh procedure can compute nonlinear crest distributions more than 25 times faster than the numerical algorithm based on the FFT technique.

In this paper we provide a perturbation analysis of finite time-inhomogeneous Markov processes. We derive closed-form representations for the derivative of the transition probability at time t, with t > 0. Elaborating on this result, we derive simple gradient estimators for transient performance characteristics either taken at some fixed point in time t, or for the integrated performance over a time interval [0 , t]. Bounds for transient performance sensitivities are presented as well. Eventually, we identify a structural property of the derivative of the generator matrix of a Markov chain that leads to a significant simplification of the estimators.

In recent years, Com–Poisson has emerged as one of the most popular discrete models in the analysis of count data owing to its flexibility in handling different types of dispersion. However, in a stationary longitudinal Com–Poisson count data set-up where the covariates are time independent, estimation of regression and dispersion parameters based on a generalized quasi-likelihood (GQL) approach involves some major computational difficulties particularly in the inversion of the joint covariance matrix. On the other hand, in practical real-life longitudinal studies, time-dependent covariates leading to non-stationary responses are more frequently encountered. This implies that further computational problems will now arise when estimating parameters under non-stationary set-ups. This paper overcomes this problem by approximating the inverse of the ill-conditioned covariance matrix in the GQL approach through a multidimensional conjugate gradient method. The performance of this novel version of the GQL approach is then assessed on simulations of AR(1) stationary and AR(1) non-stationary longitudinal Com–Poisson counts and on real-life epileptic seizure counts. However, there is not yet an algorithm to generate non-stationary longitudinal Com–Poisson counts nor a GQL algorithm to estimate the parameters under non-stationary set-ups. Thus, the paper also provides a framework to generate non-stationary AR(1) Com–Poisson counts along with the construction of a GQL equation under non-stationary set-ups.

We discuss simulation of sensitivities or Greeks of multi-asset European style options under a special Lévy process model: that is, the subordinated Brownian motion model. The Malliavin calculus method combined with Monte Carlo and quasi-Monte Carlo methods is used in the simulations. Greeks are expressed in terms of the expectations of the option payoff functions multiplied by the weights involving Malliavin derivatives for multi-asset options. Numerical results show that the Malliavin calculus method is usually more efficient than the finite difference method for options with nonsmooth payoffs. The superiority of the former method over the latter is even more significant when both are combined with quasi-Monte Carlo methods.

Suppose that X1, X2, . . . are independent identically distributed Bernoulli random variables with mean p. A Bernoulli factory for a function f takes as input X1, X2, . . . and outputs a random variable that is Bernoulli with mean f(p). A fast algorithm is a function that only depends on the values of X1, . . ., XT, where T is a stopping time with small mean. When f(p) is a real analytic function the problem reduces to being able to draw from linear functions Cp for a constant C > 1. Also it is necessary that Cp ⩽ 1 − ε for known ε > 0. Previous methods for this problem required extensive modification of the algorithm for every value of C and ε. These methods did not have explicit bounds on $\mathbb{E}[T]$ as a function of C and ε. This paper presents the first Bernoulli factory for f(p) = Cp with bounds on $\mathbb{E}[T]$ as a function of the input parameters. In fact, supp∈[0,(1−ε)/C]$\mathbb{E}[T]$ ≤ 9.5ε−1C. In addition, this method is very simple to implement. Furthermore, a lower bound on the average running time of any Cp Bernoulli factory is shown. For ε ⩽ 1/2, supp∈[0,(1−ε)/C]$\mathbb{E}[T]$≥0.004Cε−1, so the new method is optimal up to a constant in the running time.

We develop a new Monte Carlo variance reduction method to estimate the expectation of two commonly encountered path-dependent functionals: first-passage times and occupation times of sets. The method is based on a recursive approximation of the first-passage time probability and expected occupation time of sets of a Lévy bridge process that relies in part on a randomisation of the time parameter. We establish this recursion for general Lévy processes and derive its explicit form for mixed-exponential jump-diffusions, a dense subclass (in the sense of weak approximation) of Lévy processes, which includes Brownian motion with drift, Kou's double-exponential model, and hyperexponential jump-diffusion models. We present a highly accurate numerical realisation and derive error estimates. By way of illustration the method is applied to the valuation of range accruals and barrier options under exponential Lévy models and Bates-type stochastic volatility models with exponential jumps. Compared with standard Monte Carlo methods, we find that the method is significantly more efficient.

The Doss-Sussmann (DS) approach is used for uniform simulation of the Cox-Ingersoll-Ross (CIR) process. The DS formalism allows us to express trajectories of the CIR process through solutions of some ordinary differential equation (ODE) depending on realizations of a Wiener process involved. By simulating the first-passage times of the increments of the Wiener process to the boundary of an interval and solving the ODE, we uniformly approximate the trajectories of the CIR process. In this respect special attention is payed to simulation of trajectories near 0. From a conceptual point of view the proposed method gives a better quality of approximation (from a pathwise point of view) than standard, even exact, simulation of the stochastic differential equation at some deterministic time grid.

In this paper we describe a perfect simulation algorithm for the stable M/G/c queue. Sigman (2011) showed how to build a dominated coupling-from-the-past algorithm for perfect simulation of the super-stable M/G/c queue operating under first-come-first-served discipline. Sigman's method used a dominating process provided by the corresponding M/G/1 queue (using Wolff's sample path monotonicity, which applies when service durations are coupled in order of initiation of service). The method exploited the fact that the workload process for the M/G/1 queue remains the same under different queueing disciplines, in particular under the processor sharing discipline, for which a dynamic reversibility property holds. We generalise Sigman's construction to the stable case by comparing the M/G/c queue to a copy run under random assignment. This allows us to produce a naïve perfect simulation algorithm based on running the dominating process back to the time it first empties. We also construct a more efficient algorithm that uses sandwiching by lower and upper processes constructed as coupled M/G/c queues started respectively from the empty state and the state of the M/G/c queue under random assignment. A careful analysis shows that appropriate ordering relationships can still be maintained, so long as service durations continue to be coupled in order of initiation of service. We summarise statistical checks of simulation output, and demonstrate that the mean run-time is finite so long as the second moment of the service duration distribution is finite.

Let f be an integrable function on an infinite measure space (S, , π). We show that if a regenerative sequence {Xn}n≥0 with canonical measure π could be generated then a consistent estimator of λ ≡ ∫Sf dπ can be produced. We further show that under appropriate second moment conditions, a confidence interval for λ can also be derived. This is illustrated with estimating countable sums and integrals with respect to absolutely continuous measures on ℝd using a simple symmetric random walk on ℤ.

Convergence analysis is presented for recently proposed multistep schemes, when applied to a special type of forward-backward stochastic differential equations (FB-SDEs) that arises in finance and stochastic control. The corresponding k-step scheme admits a k-order convergence rate in time, when the exact solution of the forward stochastic differential equation (SDE) is given. Our analysis assumes that the terminal conditions and the FBSDE coefficients are sufficiently regular.

Electromagnetic wave propagation in complex dispersive media is governed by the time dependent Maxwell's equations coupled to equations that describe the evolution of the induced macroscopic polarization. We consider “polydispersive” materials represented by distributions of dielectric parameters in a polarization model. The work focuses on a novel computational framework for such problems involving Polynomial Chaos Expansions as a method to improve the modeling accuracy of the Debye model and allow for easy simulation using the Finite Difference Time Domain (FDTD) method. Stability and dispersion analyzes are performed for the approach utilizing the second order Yee scheme in two spatial dimensions.

In this paper we consider a classic problem concerning the high excursion probabilities of a Gaussian random field f living on a compact set T. We develop efficient computational methods for the tail probabilities ℙ{supTf(t) > b}. For each positive ε, we present Monte Carlo algorithms that run in constant time and compute the probabilities with relative error ε for arbitrarily large b. The efficiency results are applicable to a large class of Hölder continuous Gaussian random fields. Besides computations, the change of measure and its analysis techniques have several theoretical and practical indications in the asymptotic analysis of Gaussian random fields.

In this paper we study the mixing time of certain adaptive Markov chain Monte Carlo (MCMC) algorithms. Under some regularity conditions, we show that the convergence rate of importance resampling MCMC algorithms, measured in terms of the total variation distance, is O(n-1). By means of an example, we establish that, in general, this algorithm does not converge at a faster rate. We also study the interacting tempering algorithm, a simplified version of the equi-energy sampler, and establish that its mixing time is of order O(n-1/2).

We present the first class of perfect sampling (also known as exact simulation) algorithms for the steady-state distribution of non-Markovian loss systems. We use a variation of dominated coupling from the past. We first simulate a stationary infinite server system backwards in time and analyze the running time in heavy traffic. In particular, we are able to simulate stationary renewal marked point processes in unbounded regions. We then use the infinite server system as an upper bound process to simulate the loss system. The running time analysis of our perfect sampling algorithm for loss systems is performed in the quality-driven (QD) and the quality-and-efficiency-driven regimes. In both cases, we show that our algorithm achieves subexponential complexity as both the number of servers and the arrival rate increase. Moreover, in the QD regime, our algorithm achieves a nearly optimal rate of complexity.

The calculation of multivariate normal probabilities is of great importance in many statistical and economic applications. In this paper we propose a spherical Monte Carlo method with both theoretical analysis and numerical simulation. We start by writing the multivariate normal probability via an inner radial integral and an outer spherical integral using the spherical transformation. For the outer spherical integral, we apply an integration rule by randomly rotating a predetermined set of well-located points. To find the desired set, we derive an upper bound for the variance of the Monte Carlo estimator and propose a set which is related to the kissing number problem in sphere packings. For the inner radial integral, we employ the idea of antithetic variates and identify certain conditions so that variance reduction is guaranteed. Extensive Monte Carlo simulations on some probabilities confirm these claims.

Collocation has become a standard tool for approximation of parameterized systems in the uncertainty quantification (UQ) community. Techniques for least-squares regularization, compressive sampling recovery, and interpolatory reconstruction are becoming standard tools used in a variety of applications. Selection of a collocation mesh is frequently a challenge, but methods that construct geometrically unstructured collocation meshes have shown great potential due to attractive theoretical properties and direct, simple generation and implementation. We investigate properties of these meshes, presenting stability and accuracy results that can be used as guides for generating stochastic collocation grids in multiple dimensions.

This paper is concerned with the solution of the optimal stopping problem associated to the value of American options driven by continuous-time Markov chains. The value-function of an American option in this setting is characterised as the unique solution (in a distributional sense) of a system of variational inequalities. Furthermore, with continuous and smooth fit principles not applicable in this discrete state-space setting, a novel explicit characterisation is provided of the optimal stopping boundary in terms of the generator of the underlying Markov chain. Subsequently, an algorithm is presented for the valuation of American options under Markov chain models. By application to a suitably chosen sequence of Markov chains, the algorithm provides an approximate valuation of an American option under a class of Markov models that includes diffusion models, exponential Lévy models, and stochastic differential equations driven by Lévy processes. Numerical experiments for a range of different models suggest that the approximation algorithm is flexible and accurate. A proof of convergence is also provided.

We consider the stochastic Allen-Cahn equation perturbed by smooth additive Gaussian noise in a spatial domain with smooth boundary in dimension d ≤ 3, and study the semidiscretization in time of the equation by an implicit Euler method. We show that the method converges pathwise with a rate O(Δtγ) for any γ < ½. We also prove that the scheme converges uniformly in the strong Lp-sense but with no rate given.