# Applied Probability Frontiers: Computational and Modeling Challenges (15w5160)

## Organizers

Jose Blanchet (Columbia University)

Shane Henderson (Cornell University)

Donald Iglehart (Stanford University)

Thomas Kurtz (University of Wisconsin)

Pierre L'Ecuyer (Universite de Montreal)

Amy Ward (University of Southern California)

Assaf Zeevi (Columbia University)

## Objectives

The main goal of this workshop is to bring together prominent researchers in the area of applied probability (as defined in the overview section) with the following objectives:

1) Discuss the state-of-the-art on computational techniques for stochastic systems motivated by a wide range of applied problems of current interest.

2) Bring researchers from different, yet potentially synergistic areas, who have developed successful and/or promising techniques, and can contribute to generating new perspectives on problems of wide interest -- namely, beyond the boundaries of their respective area/s.

3) Discuss model features that have traditionally defined methodological advances in computation, and contrast / blend those with structures arising in new problems that stem from novel application domains

4) Identify key challenges in computational probability moving forward.

We also want to take the opportunity to celebrate the contributions to applied probability of Peter W. Glynn, whose work has significantly influenced -- and continues to influence -- key areas which are the focal points of this workshop.

### Computation of steady-state performance measures

Simulation methodology is the primary computational tool for steady-state analysis of stochastic systems. (In this document we use the terms stochastic simulation and Monte Carlo interchangeably.) There are many documented applications of the Monte Carlo method for steady-state computations, motivated by different scientific areas such as Computer Science, Operations Research, Physics, and Statistics, see [4], Ch. 4. A sampling of the main challenges here includes: acceleration of simulation schemes for exploration of a high dimensional, multi-modal, target steady state distribution, a problem that often arises in the fields of Physics and Statistics; quantification of rates of convergence to stationarity of specific Markov chains as often studied in Statistics and Computer Science; and the detection and deletion of the initial transient bias in the steady-state simulation of semi-Markov models, for which typically very little is known about the equilibrium distribution -- a problem that arises frequently in Operations Research models.

These scientific areas have witnessed exciting developments in recent years. We plan to include in the workshop a significant representation of some of these developments which, we believe can provide a basis for fruitful cross fertilization. For example, recently there has been significant progress on adaptive algorithms in the context of statistical applications (see for example [33],[26]). These algorithms are closely aligned with the theory of stochastic approximation, which is highly relevant in optimization applications arising in Operations Research. Likewise, there has been recent progress using tools traditionally studied in Operations Research, such as rare-event simulation and large deviations, to develop fast acceleration techniques for Markov chain Monte Carlo [19]. Again, these ideas are highly relevant for both communities.

These are only a couple of examples of the potential advances and cross fertilization possibilities. For example, recent advances in exact simulation (i.e. simulation without bias) open up the possibility of completely deleting the initial transient both in engineering applications and in statistics applications, which often involve unbounded and continuous state spaces [3], [12]. Finally we mentioned advances in data assimilation / filtering techniques [16], [24], which are likely to be important tools in data driven optimization settings in applied Operations Research.

### Numerical methods for SDEs and Related Continuous Time Continuous Space Models

Continuous-time models, such as stochastic (partial or even ordinary) differential equations (driven by Brownian motion or other type of noise such as Levy or long range dependent processes) are ubiquitous in scientific applications primarily because, owing to the theory of analysis (both classical and stochastic), these models are amenable to mathematical manipulation. Nevertheless, computational tasks associated with these models are often challenging given the fact that they ideally require a continuum of information to be stored in a computer's memory. Moreover, basic probabilistic quantities, such as for example, transition probability functions, are virtually impossible to compute analytically in almost all cases.

In recent years, there have been a number of major advances in the theory of simulation that have enhanced our ability to access the solution to SDEs with an increasingly accurate assessment of the control error. For example, the work of [6] has allowed us to simulate, essentially without any bias, a range of diffusion processes. This work has enabled the possibility of providing estimates for the underlying transition density and, consequently, the direct application of explicit likelihood ratio methods in statistical analysis of diffusion models [7]. Another important recent advance in numerical methods for SDEs corresponds to multilevel Monte Carlo methods, (see [20]). This approach has enabled the estimation of a large class of expectations with an optimal rate of convergence. More recently, there has also been work that allows the simulation of piecewise linear processes with sample path approximations, that are subject to a prescribed (deterministic) error with probability one [8], [9]. These constructions enable the estimation of sample path expectations of complex objects such as multidimensional local-time-like processes, for example those arising in the solution of the Skorokhod problem studied in stochastic networks [17], [9].

There are still many outstanding challenges in analyzing multidimensional diffusions, in SPDEs, and in SDEs driven by multidimensional Levy and long-range dependent processes. However, there are also opportunities for cross fertilization among the ideas discussed in the previous paragraph which have so far been studied and applied in allied areas whose communities tend to have little intersection, e.g., statistical inference (exact sampling of diffusions, [22]), mathematical finance, and chemical reaction networks (multilevel Monte Carlo, [21], [2]), and stochastic queueing networks (strong couplings, [9]).

It is important to note that a significant portion of the Markov chain Monte Carlo modeling developments discussed in the previous subsection is in continuous time [19], [26]. Thus, the workshop, having strong participation from both the steady-state Monte Carlo community and the stochastic simulation of SDEs community, provides yet another opportunity of cross fertilization on this front.

### Rare event Monte Carlo simulation

This area has been traditionally investigated in the context of communication networks, finance, and insurance. The applications related to finance and insurance have gained substantial attention in view of the recent global economic slowdown. These methods have traditionally relied on the classical theory of large deviations for light-tailed stochastic processes (see [18], [25]).

More recently, in part motivated by some of the applications described earlier, there have been advances in the development of techniques for systems with heavy-tailed characteristics [5], [25] Some of these methods rely on the application of statistical extreme value theory combined with the theory of stability of Markov chains and Lyapunov inequalities, typically used in applications involving Markov chain Monte Carlo, see [11]. The area of rare event analysis via Monte Carlo has also exported some of this techniques to the analysis of discrete structures [10] for counting applications (typically in the context of statistical analysis and computer science), and also in the analysis of Gaussian random fields and the solutions of random PDEs [1]. Moreover, some of the techniques in the area have now been applied to the design of steady-state simulation estimators [12], [9].

An important limitation in this area is that the majority of the Monte Carlo estimators that are derived and that are shown to achieve desirable optimality properties (in terms of achieving an optimal convergence rate) require the existence of asymptotic approximations. A challenge in the area is to design estimators that are both efficient in some sense and that are applicable in environments that do not require these types of approximations. In addition, it remains to be seen how the theory and algorithms developed in the rare event simulation literature can aid the development of statistical inference procedures for the likelihood of rare events, for example, using a Bayesian perspective. This, we believe, is a very promising area of potential cross fertilization between communities in Operations Research (which traditionally have investigated rare event simulation algorithms) and Statistics (who, naturally, study inference procedures).

### Simulation Optimization

Simulation optimization entails the minimization of a function that can only be evaluated through stochastic simulation. Example applications include the design of emergency-service response systems [28] and service-system staffing [14]. Algorithms for this kind of problem must cope with noisy estimates of function values and derivatives. There is a substantial literature on methods for attacking such problems; see, e.g., [4], [30] , Chapters 17-21 of [23]. The area is advancing rapidly, partly due to advances in hardware, and partly through advances in software enabled by recent research. In addition, the trend towards parallel computing, as evidenced through multicore architectures, the use of graphical processing units in scientific computing, cloud computing, and the increasing availability of high-performance computing (supercomputing), provides a rich environment for the development of new and exciting work in this area.

There are strong ties between simulation optimization and both steady-state simulation and rare-event estimation. For steady-state simulation, the method of simulated annealing with a fixed temperature generates a Markov chain, the steady-state of which is a Boltzmann distribution that concentrates near the optimal solutions of the optimization problem, e.g., [13]. The cross-entropy method exploits a link between rare-event estimation and simulation optimization, in that if one generates potential solutions from some family of distributions, then the event that one of those solutions is close to optimal is rare [32]. Furthermore, the method of stochastic approximation has already been mentioned several times in connection with the other areas of focus for the workshop, and stochastic approximation is an algorithmic device of fundamental importance in simulation optimization. A particular outstanding challenge is the design of adaptive procedures that adjust tuning sequences that calibrate the stochastic approximation algorithm, to various information that is a prior unknown concerning the target function that is to be optimized, and the environment in which it is observed/simulated (see, e.g., [15]). These adaptive techniques have a strong connection with the theory of stochastic stability in Markov chains, analysis of Markov chain Monte Carlo algorithms, and stochastic recursions (discrete time analogues of SDE's). Given these, and other links stated above, we believe there is ample potential for cross fertilization and synergies to be realized in this area as well.

Some other potential focal points include the revolution in simulation optimization algorithms made possible through parallel computing, the potential for the development of simulation optimization algorithms for enhancing methods in steady-state simulation (such as optimizing over the choice of randomized time horizon distribution in [29], [31] and rare-event estimation (such as optimizing over the choice of sampling distribution as in adaptive rare-event estimation methods).

#### Bibliography

[1] R. Adler, J. Blanchet and J. Liu (2012), “Efficient Monte Carlo for high excursions of Gaussian random field",

*Annals of Applied Probability*Vol. 22, No. 3, pp 1167-1214.

[2] D. Anderson and D. J. Higham (2012), Multi-level Monte Carlo for continuous time Markov chains, with applications in biochemical kinetics,

*SIAM: Multiscale Modeling and Simulation*Vol. 10, No. 1, 146 - 179.

[3] S. Asmussen, P. Glynn and H. Thorisson (1992), “Stationary detection in the initial transient problem”,

*ACM TOMACS*Vol 2, No. 2, pp 130-157.

[4] S. Asmussen and P. Glynn (2007),

*Stochastic Simulation: Algorithms and Analysis*1st ed. Springer, New York, NY.

[5] S. Asmussen and D. Kortscha (2013), “Error rates and improved algorithms for rare event simulation with heavy Weibull tails”,

*Methodology and Computing in Applied Probability*published online.

[6] A. Beskos and G. Roberts (2005), “Exact simulation of diffusions",

*Annals of Applied Probability*Vol. 15, No. 4, pp 2422-2444.

[7] A. Beskos, O. Papaspiliopoulos and G. Roberts (2006), “Retrospective exact simulation of diffusion sample paths with applications",

*Bernoulli*Vol. 12, No. 6, pp 1077-1098.

[8] A. Beskos, S. Peluchetti and G. Roberts (2012), “Strong simulation of the Brownian path”,

*Bernoulli*Vol. 18, No. 4, pp 1223-1248.

[9] J. Blanchet and X. Chen (2012), “Steady-state simulation of reflected Brownian motion and related stochastic networks”, http://arxiv.org/abs/1202.2062

[10] J. Blanchet (2009), “Importance Sampling and Efficient Counting for Binary Contingency Tables”,

*Annals of Applied Probability*19, pp 949 - 982.

[11] J. Blanchet and P. Glynn (2008), “Efficient Rare Event Simulation for the Single Server Queue with Heavy Tailed Distributions”,

*Annals of Applied Probability*18, 4, pp. 1351 -- 1378.

[12] J. Blanchet and K. Sigman (2011), “Perfect Sampling of Perpetuities”,

*Journal of Applied Probability*Special Vol. 48A, 165-183.

[13] P. Br'emaud (1999),

*Markov Chains: Gibbs Fields, Monte Carlo Simulation and Queues*Spinger, New York, NY.

[14] M. T. Cezik and P. L'Ecuyer (2008), “Staffing multiskill call centers via linear programming and simulation",

*Management Science*Vol. 54, No. 2, pp. 310-323.

[15] D. Cicek, M. Broadie and A. Zeevi (2012), “Multidimensional stochastic approximation: adaptive algorithms and applications",

*ACM TOMACS*to appear.

[16] Del Moral, P. (2004), “Feynman-Kac Formulae: Genealogical And Interacting Particle Systems With Applications”. Springer, New York, NY.

[17] P. Dupuis and K. Ramanan (2002), “A time-reversed representation for tail probabilities of stationary reflected Brownian motion",

*Stoch. Proc. and Their Appl.*, Vol 98, pp 253-287.

[18] P. Dupuis and H. Wang (2005), “Dynamic importance sampling for uniformly recurrent Markov Chain",

*Annals of Applied Probability*Vol. 15, pp 1-38..

[19] P. Dupuis, Y. Liu, N. Plattner and J. Doll (2012) “On the Infinite Swapping Limit for Parallel Tempering”. http://arxiv.org/abs/1110.4984.

[20] M.B. Giles (2008), “Multilevel Monte Carlo path simulation",

*Operations Research*Vol. 56, No. 3, pp 607-617.

[21] M.B. Giles and L. Szpruch (2013), “Multilevel Monte Carlo methods for applications in finance", working paper.

[22] K. Giesecke and D. Smelov (2013), “Exact sampling of jump diffusions",

*Operations Research*Vol. 61, No. 4, pp 894-907.

[23] S. G. Henderson and B. L. Nelson ed. (2006),

*Handbooks in Operations Research and Management Science*Elsevier, Amsterdam.

[24] M.A. Iglesias, K.J.H. Law and A.M. Stuart (2013), “Evaluation of Gaussian approximations for data assimilation in reservoir model",

*Computational Geoscience*Vol. 17, pp 851-885.

[25] S. Juneja (2007), “Estimating tail probabilities of heavy tailed distributions with asymptotically zero relative error",

*Queueing Systems*Vol. 57, Iss. 2-3, pp 115-127.

[26] K. Latuszynski, G.O. Roberts and J.S. Rosenthal (2013), “Adaptive Gibbs samplers and related MCMC methods",

*Annals of Applied Probability*Vol. 23, No. 1, pp 66-98.

[27] A. Mandelbaum and K. Ramanan (2010), “Directional derivatives of oblique reflection map",

*Mathematics of Operations Research*Vol. 35, pp 527-558.

[28] M. S. Maxwell and M. Restrepo and S. G. Henderson and H. Topaloglu (2010), “Approximate dynamic programming for ambulance redeployment",

*INFORMS Journal on Computing*Vol. 22, pp. 266-281.

[29] D. McLeish (2011), “A general method for debiasing a Monte Carlo estimator",

*Monte Carlo Methods and Application*Vol. 17, pp. 301-315.

[30] G. Ch. Pflug (1996),

*Optimization of Stochastic Models: The Interface Between Simulation and Optimization*Kluwer, Boston, MA.

[31] C. Rhee and P. Glynn (2012), “A new approach to unbiased estimation for SDE's", in

*Proceedings of the 2012 Winter Simulation Conference,*201-207

[32] R. Rubinstein (1999), “The cross-entropy method for combinatorial and continuous optimization",

*Methodology and Computing in Applied Probability*Vol. 1, Iss. 2, pp 127-190.

[33] A. Schreck, G. Fort and E. Moulines (2013), “Adaptive equi-energy sampler: Convergence and Illustration", working paper.