Our group studies Monte Carlo methods computational methods for random systems. The common theme in these projects is the use of mathematical analysis to understand and design new Monte Carlo methods.

- Peter Glynn, Stanford University
- Jonathan Goodman, Courant Institute, NYU
- Sangmin Lee, Courant Institute, NYU
- Jose Antonio Perez, BNP Paribas

A stochastic differential equation may be written as

A time stepping method is a way to generate sample paths that can be used for Monte Carlo. The best known methods are the (forward) Euler method and Milstein's method. Milstein's method is more accurate than the Euler method in a strong pathwise sense but not in a weaker statistical sense more appropriate for Monte Carlo. We have been exploring several related issues -- finding alternative error measures that are statistical but still measure path properties, and finding alternative time stepping methods that have the accuracy of Milstein's method but not some of its drawbacks.

In particular, we found a two stage Runge Kutta method that achieves the
accuracy of Milstein's method without using analytical derivatives of the coefficients
and without needing to sample the Levi area integral (a random variable defined by
a stochastic integral of two independent Brownian motions).
This Runge Kutta method would not be possible without our new error measure.
We are currently working on a different error measure that is closer to the
traditional strong error.
This is the expected difference between two paths, in an optimal coupling.
The so called ``coupling distance'' is related to the Wasserstein, or Monge-
Kantorovich-Rubinstein metric.

- Kranthi Gade, Merril Lynch
- Jonathan Goodman, Courant Institute, NYU

Many complex stochastic dynamical systems have subtle and very slow relaxation modes. Dynamical glasses and the glass transition are one important class of examples. We are developing a method that more accurately computes the second eigenvalue of the operator, which governs the rate at which the Markov chain statistics converges to equilibrium. So far our method is much more accurate and less noisy than Prony's method, the standard approach for these problems.

For technical details see: K. Gade, * Numerical Estimation of the Second Largest
Eigenvalue of Reversible Markov Transition Matrix *, PhD Thesis, NYU, January, 2009.

- Jonathan Goodman, Courant Institute, NYU
- Kevin Lin, University of Arizona

Control variates are a way to reduce the variance in a Monte Carlo computation by subtracting out the known answer to a simpler problem. We were interested in model one dimensional lattice systems from statistical mechanics that would satisfy detailed balance in equilibrium. Out of equilibrium, there is net energy flux through the chain that we want to estimate by Monte Carlo simulation. The problem is that long chains with a fixed temperature difference at the ends locally are close to equilibrium. This makes is hard to estimate the relatively small energy flux in the presence of statistical noise.

Our method attempts to reduce the variance by coupling the dynamic simulation to another simulation that preserves an approximation to the local equilibrium distribution. The coupling is achieved using a Metropolis style rejection strategy. Computational experiments show improvements in accuracy by roughly a factor of two, which is less than should be possible with such couplings. We hope that a more sophisticated coupling to a more accurate approximate distribution will give much better accuracy.

For technical details see: J. Goodman and K. Lin, * Dynamic control
variates for near equilibrium systems*, preprint.

- Jonathan Goodman, Courant Institute, NYU
- Sungjae Jun, Courant Institute, NYU

Sometimes it is necessary to estimate the probability of very unlikely events. For example, we may want to know the failure probability of a very reliable system. Chemical reactions are rare on the time scale of atomic motion, but do happen on a macroscopic time scale. In the case of a reliable system or a slow chemical reaction, direct simulation is not an effective way to find the rare events because most of the simulation time does not see them.

*Rare event simulation* refers to Monte Carlo methods that find rare events
more efficiently than direct simulation.
Systematic methods for rare event simulation often are based on the mathematical
theory of large deviations (for which Courant Institute mathematician S.R.S Varadhan
received the Abel Prize). Large deviation theory essentially asks the user to
find the most likely mechanism for producing rare events analytically before starting
the computation.

Our group is seeking ways to find rare events without solving the large deviation
problem first.
The picture below was produced using a histogram *bisection* strategy.
Suppose X represents a configuration of some complex random system and we
want to know the probability that f(X) > A, when this event is very small.
Our method proceeds by repeated bifurcation of the histogram. The first
step finds the median of the histogram, which is the number M1 so that the
probability of f(X)>M1 is equal to .5. The next stage is to refine this
top half to find M2 so that the probability that f(X)>m2 is .25. We do
this not by direct simulation, but by a Metropolis type sampler that rejects
configurations X that fall outside f(X)>.5. In the figure below, this
bifurcation process was repeated 20 times, giving paths whose probability
is less than one in a million. The whole process took well under a million
path simulations. The large deviation theory for this problem is known
and the path below is typical of what large deviations predicts such paths
should look like.