 
Summary: MULTILEVEL MONTE CARLO FOR CONTINUOUS TIME MARKOV CHAINS,
WITH APPLICATIONS IN BIOCHEMICAL KINETICS
DAVID F. ANDERSON1 AND DESMOND J. HIGHAM2
Abstract. We show how to extend a recently proposed multilevel Monte Carlo approach to the continuous time Markov chain setting, thereby
greatly lowering the computational complexity needed to compute expected values of functions of the state of the system to a specified accuracy.
The extension is nontrivial, exploiting a coupling of the requisite processes that is easy to simulate while providing a small variance for the
estimator. Further, and in a stark departure from other implementations of multilevel Monte Carlo, we show how to produce an unbiased estimator
that is significantly less computationally expensive than the usual unbiased estimator arising from exact algorithms in conjunction with crude Monte
Carlo. We thereby dramatically improve, in a quantifiable manner, the basic computational complexity of current approaches that have many names
and variants across the scientific literature, including the BortzKalosLebowitz algorithm, discrete event simulation, dynamic Monte Carlo, kinetic
Monte Carlo, the nfold way, the next reaction method, the residencetime algorithm, the stochastic simulation algorithm, Gillespie's algorithm,
and tauleaping. The new algorithm applies generically, but we also give an example where the coupling idea alone, even without a multilevel
discretization, can be used to improve efficiency by exploiting system structure. Stochastically modeled chemical reaction networks provide a very
important application for this work. Hence, we use this context for our notation, terminology, natural scalings, and computational examples.
Keywords: continuous time Markov chain, reaction network, computational complexity, Gillespie, next reaction
method, random time change, tauleaping, variance.
1. Introduction. This paper concerns the efficient computation of expectations for continuous time Markov
chains. Specifically, we extend the multilevel Monte Carlo approach of Giles [18], with related earlier work by
Heinrich [24], to this setting. We study the wide class of systems that can be written using the random time change
representation of Kurtz [15, Chapter 6] [32] in the form
