Monte Carlo and Quasi-Monte Carlo Methods 2008 Pierre L’Ecuyer r Art B. Owen Editors Monte Carlo and Quasi- Monte Carlo Methods 2008 Editors Pierre L’Ecuyer DIRO Université de Montréal C.P. 6128, Succ. Centre-Ville Montreal, H3C 3J7 Canada

[email protected] Art B. Owen Department of Statistics Stanford University Sequoia Hall Stanford, CA 94305 USA

[email protected] ISBN 978-3-642-04106-8 DOI 10.1007/978-3-642-04107-5 e-ISBN978-3-642-04107-5 Springer Heidelberg Dordrecht London New York Library of Congress Control Number: 2009940794 Mathematics Subject Classifi cation (2000): Primary 11K45, 65-06, 65C05, 65C10; Secondary 11K38, 65D18, 65D30, 65D32, 65R20, 91B28 ? Springer-Verlag Berlin Heidelberg 2009 This work is subject to copyright. All rights are reserved, whether the whole or part of the material is concerned, specifi cally the rights of translation, reprinting, reuse of illustrations, recitation, broadcasting, reproduction on microfi lm or in any other way, and storage in data banks. Duplication of this publication or parts thereof is permitted only under the provisions of the German Copyright Law of September 9, 1965, in its current version, and permission for use must always be obtained from Springer. Violations are liable to prosecution under the German Copyright Law. The use of general descriptive names, registered names, trademarks, etc. in this publication does not imply, evenintheabsenceofaspecifi cstatement,thatsuchnamesareexemptfromtherelevantprotective laws and regulations and therefore free for general use. Cover design: VTeX, Vilnius Printed on acid-free paper Springer is part of Springer Science+Business Media (www.springer.com) Preface This volume represents the refereed proceedings of the Eighth International Con- ference on Monte Carlo and Quasi-Monte Carlo Methods in Scientifi c Computing, which was held at the University of Montréal, from 6–11 July, 2008. It contains a limited selection of articles based on presentations made at the conference. The program was arranged with the help of an international committee consisting of: Ronald Cools, Katholieke Universiteit Leuven Luc Devroye, McGill University Henri Faure, CNRS Marseille Paul Glasserman, Columbia University Peter W. Glynn, Stanford University Stefan Heinrich, University of Kaiserslautern Fred J. Hickernell, Illinois Institute of Technology Aneta Karaivanova, Bulgarian Academy of Science Alexander Keller, mental images GmbH, Berlin Adam Kolkiewicz, University of Waterloo Frances Y. Kuo, University of New South Wales Christian Lécot, Université de Savoie, Chambéry Pierre L’Ecuyer, Université de Montréal (Chair and organizer) Jun Liu, Harvard University Peter Mathé, Weierstrass Institute Berlin Makoto Matsumoto, Hiroshima University Thomas Müller-Gronbach, Otto von Guericke Universit?t Harald Niederreiter, National University of Singapore Art B. Owen, Stanford University Gilles Pagès, Université Pierre et Marie Curie (Paris 6) Klaus Ritter, TU Darmstadt Karl Sabelfeld, Weierstrass Institute Berlin Wolfgang Ch. Schmid, University of Salzburg Ian H. Sloan, University of New South Wales Jerome Spanier, University of California, Irvine Bruno Tuffi n, IRISA-INRIA, Rennes Henryk Wo′ zniakowski, Columbia University. v viPreface The local arrangements (program production, publicity, web site, registration, social events, etc.) were ably handled by Carole Dufour (GERAD), Marilyne Lavoie (GERAD), Louis Pelletier (CRM), Marie Perreault (GERAD), and Suzette Paradis (CRM). Francine Benoit (GERAD) helped with editing the proceedings. This conference continued the tradition of biennial MCQMC conferences initi- ated by Harald Niederreiter. They were begun at the University of Nevada in Las Vegas, Nevada, USA, in June 1994 and followed by conferences at the University of Salzburg, Austria, in July 1996, the Claremont Colleges in Claremont, Califor- nia, USA, in June 1998, Hong Kong Baptist University in Hong Kong, China, in November 2000, the National University of Singapore, Republic of Singapore, in November 2002, the Palais des Congrès in Juan-les-Pins, France, in June 2004, and Ulm University, Germany, in July 2006. The next MCQMC conference will be held in Warsaw, Poland, in August 2010. The proceedings of these previous conferences were all published by Springer- Verlag, under the following titles: ? MonteCarloandQuasi-MonteCarloMethodsinScientifi cComputing(H.Nieder- reiter and P.J.-S. Shiue, eds.), ?Monte Carlo and Quasi-Monte Carlo Methods 1996 (H. Niederreiter, P. Hellekalek, G. Larcher and P. Zinterhof, eds.), ?Monte Carlo and Quasi-Monte Carlo Methods 1998 (H. Niederreiter and J. Spanier, eds.), ?Monte Carlo and Quasi-Monte Carlo Methods 2000 (K.-T. Fang, F.J. Hickernell and H. Niederreiter, eds.), ?Monte Carlo and Quasi-Monte Carlo Methods 2002 (H. Niederreiter, ed.), ?Monte Carlo and Quasi-Monte Carlo Methods 2004 (H. Niederreiter and D. Ta- lay, eds.), ?Monte Carlo and Quasi-Monte Carlo Methods 2006 (A. Keller and S. Heinrich and H. Niederreiter, eds.). The program of the conference was rich and varied with over 135 talks being presented. Highlights were the invited plenary talks given by Josef Dick (Univer- sity of New South Wales), Arnaud Doucet (University of British Columbia), Daan Frenkel (University of Cambridge), Paul Glasserman (Columbia University), Chris- tiane Lemieux (University of Waterloo), Jun Liu (Harvard University), Klaus Ritter (TU Darmstadt), Jeffrey Rosenthal (University of Toronto), Wolfgang Schmid (Uni- versity of Salzburg), and Andrew Stuart (Warwick University). The papers in this volume were carefully screened and cover both the theory and the applications of Monte Carlo and quasi-Monte Carlo methods. We thank the anonymous reviewers for their reports and many others who con- tributed enormously to the excellent quality of the conference presentations and to the high standards for publication in these proceedings by careful review of the ab- stracts and manuscripts that were submitted. We gratefully acknowledge generous fi nancial support of the conference by the Centre de Recherches Mathématiques (CRM), the Groupe d’études et de Recherche en Analyse de Décisions (GERAD), Mathematics for Information Technology Prefacevii and Complex Systems (MITACS), and the American National Science Foundation (NSF). Finally, we want to express our gratitude to Springer-Verlag for publishing this volume. Pierre L’Ecuyer July 2009Art Owen Contents Part I Tutorials Monte Carlo and Quasi-Monte Carlo for Statistics .3 Art B. Owen Monte Carlo Computation in Finance 19 Jeremy Staum Part II Invited Articles Particle Markov Chain Monte Carlo for Effi cient Numerical Simulation .45 Christophe Andrieu, Arnaud Doucet, and Roman Holenstein Computational Complexity of Metropolis-Hastings Methods in High Dimensions.61 Alexandros Beskos and Andrew Stuart On Quasi-Monte Carlo Rules Achieving Higher Order Convergence 73 Josef Dick Sensitivity Estimates for Compound Sums 97 Paul Glasserman and Kyoung-Kuk Kim New Perspectives on (0,s)-Sequences . 113 Christiane Lemieux and Henri Faure Variable Subspace Sampling and Multi-level Algorithms 131 Thomas Müller-Gronbach and Klaus Ritter Markov Chain Monte Carlo Algorithms: Theory and Practice . 157 Jeffrey S. Rosenthal ix xContents MINT – New Features and New Results. 171 Rudolf Schürer and Wolfgang Ch. Schmid Part III Contributed Articles Recursive Computation of Value-at-Risk and Conditional Value-at-Risk using MC and QMC . 193 Olivier Bardou, Noufel Frikha, and Gilles Pagès Adaptive Monte Carlo Algorithms Applied to Heterogeneous Transport Problems 209 Katherine Bhan, Rong Kong, and Jerome Spanier Effi cient Simulation of Light-Tailed Sums: an Old-Folk Song Sung to a Faster New Tune. 227 Jose H. Blanchet, Kevin Leder, and Peter W. Glynn Distribution of Digital Explicit Inversive Pseudorandom Numbers and Their Binary Threshold Sequence . 249 Zhixiong Chen, Domingo Gomez, and Arne Winterhof Extensions of Fibonacci Lattice Rules 259 Ronald Cools and Dirk Nuyens Effi cient Search for Two-Dimensional Rank-1 Lattices with Applications in Graphics 271 Sabrina Dammertz, Holger Dammertz, and Alexander Keller Parallel Random Number Generators Based on Large Order Multiple Recursive Generators 289 Lih-Yuan Deng, Jyh-Jen Horng Shiau, and Gwei-Hung Tsai Effi cient Numerical Inversion for Financial Simulations . 297 Gerhard Derfl inger, Wolfgang H?rmann, Josef Leydold, and Halis Sak Equidistribution Properties of Generalized Nets and Sequences 305 Josef Dick and Jan Baldeaux Implementation of a Component-By-Component Algorithm to Generate Small Low-Discrepancy Samples 323 Benjamin Doerr, Michael Gnewuch, and Magnus Wahlstr?m Quasi-Monte Carlo Simulation of Diffusion in a Spatially Nonhomogeneous Medium 339 Rami El Haddad, Christian Lécot, and Gopalakrishnan Venkiteswaran L2Discrepancy of Two-Dimensional Digitally Shifted Hammersley Point Sets in Base b . 355 Henri Faure and Friedrich Pillichshammer Contentsxi Vibrato Monte Carlo Sensitivities 369 Michael B. Giles The Weighted Variance Minimization in Jump-Diffusion Stochastic Volatility Models 383 Anatoly Gormin and Yuri Kashtanov (t,m,s)-Nets and Maximized Minimum Distance, Part II 395 Leonhard Grünschlo? and Alexander Keller Automation of Statistical Tests on Randomness to Obtain Clearer Conclusion . 411 Hiroshi Haramoto On Subsequences of Niederreiter-Halton Sequences 423 Roswitha Hofer Correcting the Bias in Monte Carlo Estimators of American-style Option Values 439 K.H. Felix Kan, R. Mark Reesor, Tyson Whitehead, and Matt Davison Fast Principal Components Analysis Method for Finance Problems With Unequal Time Steps 455 Jens Keiner and Benjamin J. Waterhouse Adaptive Monte Carlo Algorithms for General Transport Problems . 467 Rong Kong, Martin Ambrose, and Jerome Spanier On Array-RQMC for Markov Chains: Mapping Alternatives and Convergence Rates 485 Pierre L’Ecuyer, Christian Lécot, and Adam L’Archevêque-Gaudet Testing the Tests: Using Random Number Generators to Improve Empirical Tests . 501 Paul Leopardi Stochastic Spectral Formulations for Elliptic Problems 513 Sylvain Maire and Etienne Tanré Adaptive (Quasi-)Monte Carlo Methods for Pricing Path-Dependent Options 529 Roman N. Makarov Monte Carlo Simulation of Stochastic Integrals when the Cost of Function Evaluation Is Dimension Dependent . 545 Ben Niu and Fred J. Hickernell xiiContents Recent Progress in Improvement of Extreme Discrepancy and Star Discrepancy of One-Dimensional Sequences . 561 Victor Ostromoukhov Discrepancy of Hyperplane Nets and Cyclic Nets . 573 Friedrich Pillichshammer and Gottlieb Pirsic A PRNG Specialized in Double Precision Floating Point Numbers Using an Affi ne Transition . 589 Mutsuo Saito and Makoto Matsumoto On the Behavior of the Weighted Star Discrepancy Bounds for Shifted Lattice Rules . 603 Vasile Sinescu and Pierre L’Ecuyer Ergodic Estimations of Upscaled Coeffi cients for Diffusion in Random Velocity Fields 617 Nicolae Suciu and C? alin Vamo? s Green’s Functions by Monte Carlo. 627 David White and Andrew Stuart Tractability of Multivariate Integration for Weighted Korobov Spaces: My 15 Year Partnership with Ian Sloan . 637 Henryk Wo′ zniakowski Conference Participants. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 655 Index . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 671 Part I Tutorials Monte Carlo and Quasi-Monte Carlo for Statistics Art B. Owen Abstract This article reports on the contents of a tutorial session at MCQMC 2008. The tutorial explored various places in statistics where Monte Carlo methods can be used. There was a special emphasis on areas where Quasi-Monte Carlo ideas have been or could be applied, as well as areas that look like they need more research. 1 Introduction This survey is aimed at exposing good problems in statistics to researchers in Quasi- Monte Carlo. It has a mix of well known and not so well known topics, which both have their place in a research context. The selection of topics is tilted in the direction of problems that I have looked at. That enables me to use real examples, and examples are crucial to understanding statistics. Monte Carlo methods are ubiquitous in statistics. Section 2 presents the boot- strap. It is a method of resampling the observed data to judge the uncertainty in a quantity. The bootstrap makes minimal assumptions about how the data were ob- tained. Some efforts at bringing balance to the resampling process have brought improvements, but they are not large enough to have made much impact on how the bootstrap is used. Permutation tests, considered in Section 3 have a similar fl avor to the bootstrap, but there, efforts to impose balance can distort the results. Markov chain Monte Carlo (Section 4) is used when we cannot directly sample the quantity of interest, but are at least able to fi nd a Markov chain from whose stationary distribution the desired quantity can be sampled. Space limitations make it impossible to cover all of the topics from a three hour tutorial in depth. The work on QMC for MCMC has appeared in [31], [43] and in Tribble’s dissertation [42], and so it is just sketched here. Department of Statistics, Stanford University, Stanford CA, 94305 url: //stat.stanford.edu/?owen P. L’Ecuyer, A.B. Owen (eds.), Monte Carlo and Quasi-Monte Carlo Methods 2008, DOI 10.1007/978-3-642-04107-5 1, ? Springer-Verlag Berlin Heidelberg 2009 3 4Art B. Owen Fig. 1 A measure of kidney damage is plotted versus age for 60 subjects in [34]. Monte Carlo methods are used for search as well as for integration. Section 5 presents the method of least trimmed squares (LTS). This is the most effective method known for highly robust regression model fi tting. It uses an ad hoc Monte Carlo search strategy. Quasi-Monte Carlo methods have been used in search prob- lems [26, Chapter 6], usually to search over the unit cube [27]. The space to search over in LTS is of a combinatorial nature. Finally, Section 6 points to two more im- portant problems from statistics where QMC may be useful. 2 The Bootstrap Figure 1 plots a measure Yiof kidney damage against age Xi, for 60 subjects studied in Rodwell et al. [34]. There is a clear tendency for older subjects to have greater kidney damage. This may be quantifi ed through the correlation coeffi cient, which on this data takes the value 0.59. Since only 60 subjects were measured, we very much doubt that the true correlation taken over all people is exactly 0.59. Using ρ to denote a hypothetical true correlation and ? ρ to denote the measured one, we may just want to know the variance V( ? ρ) of our estimate. The variance of the sample correlation is not hard to fi nd in closed form, so long as the data are a sample from the bivariate normal distribution. But we have no reason to expect that assumption is good enough to use. The bootstrap, introduced by Efron [5] provides a way out of that assumption. Operationally one does the following: 1. For b = 1,.,B MC and QMC for Statistics5 2. Draw (X?b i ,Y?b i ), 1 ≤ i ≤ 60 with replacement from the original data. 3. Compute ? ρ?b= ? ρ?(X?b 1 ,Y?b 1 ),.,(X?b 60,Y ? 60) ?. 4. Return the variance of ? ρ?1,. ? ρ?B. Using B = 9999 the result came out to be 0.0081, so that the standard deviation of the ? ρ?values is 0.090. What we actually got was a Monte Carlo estimate of the variance of the bootstrapped correlations ? ρ?when resampling from the data. Even if we let B → ∞ this would not be the same as the variance we want, which is that of ? ρ when sampling from the unknown true distribution of (X,Y) pairs. But bootstrap theory shows that the two variances become close quickly as the number n of sample values increases [7]. First impressions of the bootstrap are either that it is obviously ok, or that it is somehow too good to be true, like pulling oneself up by the bootstraps. The formal justifi cation of the bootstrap begins with a statistical functional T defi ned on dis- tributions F. In this case T(F) gives the variance of the correlation measured on 60 pairs of data drawn from the distribution F on R2. Let F0be the unknown true distribution and?Fnbe the distribution that puts equal probability 1/n on all n data points. As n increases?Fnapproaches F0. Then a continuity argument gives T(?Fn) approaching T(F0). The continuity argument holds in great generality but there are exceptions, as well as remedies in some of those cases [32]. The bootstrap can also be used to estimate and correct for biases in statistics. Let E( ? ρ | F) denote the expected value of ? ρ when sampling n data pairs from F. Typically E( ? ρ) ? = ρ, so that the sample correlation is biased. The bootstrap estimate ofthebiasB(F)≡E( ? ρ |F)?ρ(F)isB(?Fn)≡E( ? ρ?|?Fn)?ρ(?Fn).Wecanestimate this bias by resampling. In the present example we fi nd that the average value of ? ρ?? ? ρ in resampling is ?0.0047. If we are worried that ? ρ is too small by 0.0047 we can add 0.0047 (i.e. subtract the estimated bias) and get 0.594 instead of 0.589. Here the bias adjustment is very small. That is typical unless the method used has a large number of parameters relative to the sample size. Figure 2 shows a histogram of the 9,999 bootstrap samples used in this anal- ysis. The histogram is skewed, centered near the original correlation, and is quite wide. The resampled correlations cut the real line into 10,000 intervals. A bootstrap 95% confi dence interval is formed by taking the central 9500 of those intervals. If the values are sorted ? ρ?(1)≤ ? ρ?(2)≤ ··· ≤ ? ρ?(9999) , then the 95% confi dence in- terval goes from ? ρ?(250)to ? ρ?(9750) . In this example we have 95% confi dence that 0.391≤ρ ≤0.755. Similarly there is 99% confi dence that ? ρ?(100)≤ρ ≤ ? ρ?(9900), or 0.346≤ρ ≤0.765. Bootstrap confi dence levels are not exact. They are approximate confi denceintervals.Typicallytheyhavecoverageprobabilityequaltotheirnominal level plus O(n?1). See [13]. The intervals presented here are known as percentile intervals. They are the simplest but not the only bootstrap confi dence interval. See [7] for other choices. The balanced bootstrap [4] is an attempt to improve on bootstrap re-sampling. Instead of sampling n observations with replacement B times, it forms a large pool of nB observations, with B copies of each original data point. Then it randomly par- titions them into B subsets of equal size n. Those groups are treated as the bootstrap 6Art B. Owen Fig. 2 This fi gure shows 9999 bootstrap resampled correlations for the kidney data. There is a reference point at the sample correlation. Solid vertical lines enclose the central 95% of the histogram. Dashed lines enclose the central 99%. samples. Now each original observation appears the same number B of times among the sampled data. The balanced bootstrap is similar to QMC. Higher order balanc- ing, based on orthogonal arrays, was proposed by [12], but that proposal requires the number n of observations to be a prime power. To apply QMC, it helps to frame the bootstrap problem as an integral over[0,1]n, as discussed in [28] and thoroughly investigated by Liu [21]. We may write X? i as X?nUi?where Ui～ U(0,1) are independent. Then X?= (X? 1,.,X ? n) is a function of U ～ U(0,1)n. The bootstrap estimate of bias is a Monte Carlo estimate on B samples of ? [0,1]nQ(U)dU for Q(U) = T(X?(U))?T(X),(1) where T(X) is a shorthand for T(n?1 ?n i=1δXi) with δx the point mass distribution at x. The bootstrap estimate of variance has integrand Q(U) = ? T(X?(U))? ? [0,1]n T(X?(U))dU ?2 .(2) The upper end of a bootstrap 95% confi dence interval is the solution T 0.975 to ? [0,1]nQ(U)dU = 0.975 where Q(U) = 1T(X?(U))≤T0.975,(3) while the lower end uses 0.025 instead of 0.975. MC and QMC for Statistics7 To implement the ordinary bootstrap we take points Ub= (Ub1,.,Ubn) ～ U(0,1)nfor b = 1,.,B and use ordinary Monte Carlo estimates of the