Multilevel Monte Carlo Research Articles

An explicit positivity-preserving scheme for the Heston 3/2-model with order-one strong convergence

This article is concerned with numerical approximations of the Heston 3/2-model from mathematical finance, which takes values in (0,∞) and possesses superlinearly growing drift and diffusion coefficients. To discretize the SDE model, a new Milstein-type scheme is proposed, which can be explicitly solved and is positivity-preserving unconditionally, i.e., for any time step-size h>0. Furthermore, a mean-square convergence rate of order one is proved in the non-globally Lipschitz regime, which is highly non-trivial, by noting that the diffusion coefficient grows super-linearly. The above theoretical results can be then used to justify the multilevel Monte Carlo (MLMC) methods for approximating expectations of some functions of the solution to the Heston 3/2-model. Indeed, the unconditional positivity-preserving property is particularly desirable in the MLMC setting, where large discretization time steps are used. The obtained order-one convergence in turn promises the desired relevant variance of the multilevel estimator and justifies the optimal complexity O(ϵ−2) for the MLMC approach, where ϵ>0 is the required target accuracy. Numerical experiments are finally reported to confirm the above results.

A Weak MLMC Scheme for Lévy-Copula-Driven SDEs with Applications to the Pricing of Credit, Equity and Interest Rate Derivatives

This paper develops a novel weak multilevel Monte-Carlo (MLMC) approximation scheme for Lévy-driven Stochastic Differential Equations (SDEs). The scheme is based on the state space discretization (via a continuous-time Markov chain approximation [Mijatović, Aleksandar, Matija Vidmar, and Saul Jacka. 2012. “Markov Chain Approximations for Transition Densities of Lévy Processes.” Electronic Journal of Probability 19]) of the pure-jump component of the driving Lévy process and is particularly suited if the multidimensional driver is given by a Lévy copula. The multilevel version of the algorithm requires a new coupling of the approximate Lévy drivers in the consecutive levels of the scheme, which is defined via a coupling of the corresponding Poisson point processes. The multilevel scheme is weak in the sense that the bound on the level variances is based on the coupling alone without requiring strong convergence. Moreover, the coupling is natural for the proposed discretization of jumps and is easy to simulate. The approximation scheme and its multilevel analogous are applied to examples taken from mathematical finance, including the pricing of credit, equity and interest rate derivatives.

Multilevel approximation of Gaussian random fields: Covariance compression, estimation, and spatial prediction

The distribution of centered Gaussian random fields (GRFs) indexed by compacta such as smooth, bounded Euclidean domains or smooth, compact and orientable manifolds is determined by their covariance operators. We consider centered GRFs given as variational solutions to coloring operator equations driven by spatial white noise, with an elliptic self-adjoint pseudodifferential coloring operator from the Hörmander class. This includes the Matérn class of GRFs as a special case. Using biorthogonal multiresolution analyses on the manifold, we prove that the precision and covariance operators, respectively, may be identified with bi-infinite matrices and finite sections may be diagonally preconditioned rendering the condition number independent of the dimension p of this section. We prove that a tapering strategy by thresholding applied on finite sections of the bi-infinite precision and covariance matrices results in optimally numerically sparse approximations. That is, asymptotically only linearly many nonzero matrix entries are sufficient to approximate the original section of the bi-infinite covariance or precision matrix using this tapering strategy to arbitrary precision. The locations of these nonzero matrix entries can be determined a priori. The tapered covariance or precision matrices may also be optimally diagonally preconditioned. Analysis of the relative size of the entries of the tapered covariance matrices motivates novel, multilevel Monte Carlo (MLMC) oracles for covariance estimation, in sample complexity that scales log-linearly with respect to the number p of parameters. In addition, we propose and analyze novel compressive algorithms for simulating and kriging of GRFs. The complexity (work and memory vs. accuracy) of these three algorithms scales near-optimally in terms of the number of parameters p of the sample-wise approximation of the GRF in Sobolev scales.

Antithetic multilevel Monte Carlo method for approximations of SDEs with non-globally Lipschitz continuous coefficients

In the field of computational finance, one is commonly interested in the expected value of a financial derivative whose payoff depends on the solution of stochastic differential equations (SDEs). For multi-dimensional SDEs with non-commutative diffusion coefficients in the globally Lipschitz setting, a kind of one-half order truncated Milstein-type scheme without Lévy areas was recently introduced by Giles and Szpruch (2014), which combined with the antithetic multilevel Monte Carlo (MLMC) gives the optimal overall computational cost O(ϵ−2) for the required target accuracy ϵ. Nevertheless, many nonlinear SDEs in applications have non-globally Lipschitz continuous coefficients and the corresponding theoretical guarantees for antithetic MLMC are absent in the literature. In the present work, we aim to fill the gap and analyze antithetic MLMC in a non-globally Lipschitz setting. First, we propose a family of modified Milstein-type schemes without Lévy areas to approximate SDEs with non-globally Lipschitz continuous coefficients. The expected one-half order of strong convergence is recovered in a non-globally Lipschitz setting, where even the diffusion coefficients are allowed to grow superlinearly. This then helps us to analyze the relevant variance of the multilevel estimator and the optimal computational cost is finally achieved for the antithetic MLMC. Since getting rid of the Lévy areas destroys the martingale properties of the scheme, the analysis of both the convergence rate and the desired variance becomes highly non-trivial in the non-globally Lipschitz setting. By introducing an auxiliary approximation process, we develop non-standard arguments to overcome the essential difficulties. Numerical experiments are provided to confirm the theoretical findings.

Robust Average-Reward Reinforcement Learning

Robust Markov decision processes (MDPs) aim to find a policy that optimizes the worst-case performance over an uncertainty set of MDPs. Existing studies mostly have focused on the robust MDPs under the discounted reward criterion, leaving the ones under the average-reward criterion largely unexplored. In this paper, we develop the first comprehensive and systematic study of robust average-reward MDPs, where the goal is to optimize the long-term average performance under the worst case. Our contributions are four-folds: (1) we prove the uniform convergence of the robust discounted value function to the robust average-reward function as the discount factor γ goes to 1; (2) we derive the robust average-reward Bellman equation, characterize the structure of its solution set, and prove the equivalence between solving the robust Bellman equation and finding the optimal robust policy; (3) we design robust dynamic programming algorithms, and theoretically characterize their convergence to the optimal policy; and (4) we design two model-free algorithms unitizing the multi-level Monte-Carlo approach, and prove their asymptotic convergence

Multilevel Monte Carlo Methods for Stochastic Convection–Diffusion Eigenvalue Problems

We develop new multilevel Monte Carlo (MLMC) methods to estimate the expectation of the smallest eigenvalue of a stochastic convection–diffusion operator with random coefficients. The MLMC method is based on a sequence of finite element (FE) discretizations of the eigenvalue problem on a hierarchy of increasingly finer meshes. For the discretized, algebraic eigenproblems we use both the Rayleigh quotient (RQ) iteration and implicitly restarted Arnoldi (IRA), providing an analysis of the cost in each case. By studying the variance on each level and adapting classical FE error bounds to the stochastic setting, we are able to bound the total error of our MLMC estimator and provide a complexity analysis. As expected, the complexity bound for our MLMC estimator is superior to plain Monte Carlo. To improve the efficiency of the MLMC further, we exploit the hierarchy of meshes and use coarser approximations as starting values for the eigensolvers on finer ones. To improve the stability of the MLMC method for convection-dominated problems, we employ two additional strategies. First, we consider the streamline upwind Petrov–Galerkin formulation of the discrete eigenvalue problem, which allows us to start the MLMC method on coarser meshes than is possible with standard FEs. Second, we apply a homotopy method to add stability to the eigensolver for each sample. Finally, we present a multilevel quasi-Monte Carlo method that replaces Monte Carlo with a quasi-Monte Carlo (QMC) rule on each level. Due to the faster convergence of QMC, this improves the overall complexity. We provide detailed numerical results comparing our different strategies to demonstrate the practical feasibility of the MLMC method in different use cases. The results support our complexity analysis and further demonstrate the superiority over plain Monte Carlo in all cases.

Uncertainty quantification in the Henry problem using the multilevel Monte Carlo method

We investigate the applicability of the well-known multilevel Monte Carlo (MLMC) method to the class of density-driven flow problems, in particular the problem of salinisation of coastal aquifers. As a test case, we solve the uncertain Henry saltwater intrusion problem. Unknown porosity, permeability and recharge parameters are modelled by using random fields. The classical deterministic Henry problem is non-linear and time-dependent, and can easily take several hours of computing time. Uncertain settings require the solution of multiple realisations of the deterministic problem, and the total computational cost increases drastically. Instead of computing of hundreds random realisations, typically the mean value and the variance are computed. The standard methods such as the Monte Carlo or surrogate-based methods are a good choice, but they compute all stochastic realisations on the same, often, very fine mesh. They also do not balance the stochastic and discretisation errors. These facts motivated us to apply the MLMC method. We demonstrate that by solving the Henry problem on multi-level spatial and temporal meshes, the MLMC method reduces the overall computational and storage costs. To reduce the computing cost further, parallelization is performed in both physical and stochastic spaces. To solve each deterministic scenario, we run the parallel multigrid solver ug4 in a black-box fashion.

Efficient risk estimation via nested multilevel quasi-Monte Carlo simulation

We consider the problem of estimating the probability of a large loss from a financial portfolio, where the future loss is expressed as a conditional expectation. Since the conditional expectation is intractable in most cases, one may resort to nested simulation. To reduce the complexity of nested simulation, we present an improved multilevel Monte Carlo (MLMC) method by using quasi-Monte Carlo (QMC) to estimate the portfolio loss in each financial scenario generated via Monte Carlo. We prove that using QMC can accelerate the convergence rates in both the crude nested simulation and the multilevel nested simulation. Under certain conditions, the complexity of the proposed MLMC method can be reduced to O(ϵ−2(logϵ)2). On the other hand, we find that using QMC in MLMC encounters a high-kurtosis phenomenon due to the existence of indicator functions. To remedy this, we propose a smoothed method which uses logistic sigmoid functions to approximate indicator functions. Numerical results show that the optimal MLMC complexity O(ϵ−2) is almost attained even in moderate high dimensions.

MULTILEVEL MONTE CARLO ESTIMATORS FOR DERIVATIVE-FREE OPTIMIZATION UNDER UNCERTAINTY

Optimization is a key tool for scientific and engineering applications; however, in the presence of models affected by uncertainty, the optimization formulation needs to be extended to consider statistics of the quantity of interest. Optimization under uncertainty (OUU) deals with this endeavor and requires uncertainty quantification analyses at several design locations; i.e., its overall computational cost is proportional to the cost of performing a forward uncertainty analysis at each design location. An OUU workflow has two main components: an inner loop strategy for the computation of statistics of the quantity of interest, and an outer loop optimization strategy tasked with finding the optimal design, given a merit function based on the inner loop statistics. In this work, we propose to alleviate the cost of the inner loop uncertainty analysis by leveraging the so-called multilevel Monte Carlo (MLMC) method, which is able to allocate resources over multiple models with varying accuracy and cost. The resource allocation problem in MLMC is formulated by minimizing the computational cost given a target variance for the estimator. We consider MLMC estimators for statistics usually employed in OUU workflows and solve the corresponding allocation problem. For the outer loop, we consider a derivative-free optimization strategy implemented in the SNOWPAC library; our novel strategy is implemented and released in the Dakota software toolkit. We discuss several numerical test cases to showcase the features and performance of our approach with respect to its Monte Carlo single fidelity counterpart.

Gradient-based optimisation of the conditional-value-at-risk using the multi-level Monte Carlo method

In this work, we tackle the problem of minimising the Conditional-Value-at-Risk (CVaR) of output quantities of complex differential models with random input data, using gradient-based approaches in combination with the Multi-Level Monte Carlo (MLMC) method. In particular, we consider the framework of multi-level Monte Carlo for parametric expectations introduced in [24] and propose modifications of the MLMC estimator, error estimation procedure, and adaptive MLMC parameter selection to ensure the estimation of the CVaR and sensitivities for a given design with a prescribed accuracy. We then propose combining the MLMC framework with an alternating inexact minimisation-gradient descent algorithm, for which we prove Q-linear convergence in the optimisation iterations under the assumptions of strong convexity and Lipschitz continuity of the gradient of the objective function. We demonstrate the performance of our approach on two numerical examples of practical relevance, which evidence the same optimal asymptotic cost-tolerance behaviour as standard MLMC methods for fixed design computations of output expectations.

Multilevel Monte Carlo FEM for elliptic PDEs with Besov random tree priors

We develop a multilevel Monte Carlo (MLMC)-FEM algorithm for linear, elliptic diffusion problems in polytopal domain D⊂Rd\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$${\\mathcal {D}}\\subset {\\mathbb {R}}^d$$\\end{document}, with Besov-tree random coefficients. This is to say that the logarithms of the diffusion coefficients are sampled from so-called Besov-tree priors, which have recently been proposed to model data for fractal phenomena in science and engineering. Numerical analysis of the fully discrete FEM for the elliptic PDE includes quadrature approximation and must account for (a) nonuniform pathwise upper and lower coefficient bounds, and for (b) low path-regularity of the Besov-tree coefficients. Admissible non-parametric random coefficients correspond to random functions exhibiting singularities on random fractals with tunable fractal dimension, but involve no a-priori specification of the fractal geometry of singular supports of sample paths. Optimal complexity and convergence rate estimates for quantities of interest and for their second moments are proved. A convergence analysis for MLMC-FEM is performed which yields choices of the algorithmic steering parameters for efficient implementation. A complexity (“error vs work”) analysis of the MLMC-FEM approximations is provided.

Uncertainty quantification in coastal aquifers using the multilevel Monte Carlo method

AbstractWe are solving a problem of salinisation of coastal aquifers. As a test case example, we consider the Henry saltwater intrusion problem. Since porosity, permeability and recharge are unknown or only known at a few points, we model them using random fields and random variables. The Henry problem describes a two‐phase flow and is non‐linear and time‐dependent. The solution to be found is the expectation of the salt mass fraction, which is uncertain and time‐dependent. To estimate this expectation, we use the well‐known multilevel Monte Carlo (MLMC) method. The MLMC method takes just a few samples on computationally expensive (fine) meshes and more samples on cheap (coarse) meshes. Then, by building a telescoping sum, the MLMC method estimates the expected value at a much lower computational cost than the classical Monte Carlo method. The deterministic solver used here is the well‐known parallel and scalable ug4 solver.

Optimizing shift selection in multilevel Monte Carlo for disconnected diagrams in lattice QCD

The calculation of disconnected diagram contributions to physical signals is a computationally expensive task in Lattice QCD. To extract the physical signal, the trace of the inverse Lattice Dirac operator, a large sparse matrix, must be stochastically estimated. Because the variance of the stochastic estimator is typically large, variance reduction techniques must be employed. Multilevel Monte Carlo (MLMC) methods reduce the variance of the trace estimator by utilizing a telescoping sequence of estimators. Frequency Splitting is one such method that uses a sequence of inverses of shifted operators to estimate the trace of the inverse lattice Dirac operator, however there is no a priori way to select the shifts that minimize the cost of the multilevel trace estimation. In this article, we present a sampling and interpolation scheme that is able to predict the variances associated with Frequency Splitting under displacements of the underlying space time lattice. The interpolation scheme is able to predict the variances to high accuracy and therefore chooses shifts that correspond to an approximate minimum of the cost for the trace estimation. We show that Frequency Splitting with the chosen shifts displays significant speedups over multigrid deflation, and that these shifts can be used for multiple configurations within the same ensemble with no penalty to performance.

A massively parallel implementation of multilevel Monte Carlo for finite element models

The Multilevel Monte Carlo (MLMC) method has proven to be an effective variance-reduction statistical method for Uncertainty Quantification (UQ) in Partial Differential Equation (PDE) models, combining model computations at different levels to create an accurate estimate. Still, the computational complexity of the resulting method is extremely high, particularly for 3D models, which requires advanced algorithms for the efficient exploitation of High Performance Computing (HPC). In this article we present a new implementation of the MLMC in massively parallel computer architectures, exploiting parallelism within and between each level of the hierarchy. The numerical approximation of the PDE is performed using the finite element method but the algorithm is quite general and could be applied to other discretization methods. The two key ingredients of the implementation are a good processor partition scheme together with a good scheduling algorithm to assign work to different processors. We introduce a multiple partition of the set of processors that permits the simultaneous execution of different levels and we develop a dynamic scheduling algorithm to exploit it. The problem of finding the optimal scheduling of distributed tasks in a parallel computer is an NP-complete problem. We propose and analyze a new greedy scheduling algorithm to assign samples and we show that it is a 2-approximation, which is the best that may be expected under general assumptions. On top of this result we design a distributed memory implementation using the Message Passing Interface (MPI) standard. Finally we present a set of numerical experiments illustrating its scalability properties.

Urban Flood Modeling: Uncertainty Quantification and Physics‐Informed Gaussian Processes Regression Forecasting

AbstractEstimating uncertainty in flood model predictions is important for many applications, including risk assessment and flood forecasting. We focus on uncertainty in physics‐based urban flooding models. We consider the effects of the model's complexity and uncertainty in key input parameters. The effect of rainfall intensity on the uncertainty in water depth predictions is also studied. As a test study, we choose the Interconnected Channel and Pond Routing (ICPR) model of a part of the city of Minneapolis. The uncertainty in the ICPR model's predictions of the floodwater depth is quantified in terms of the ensemble variance using the multilevel Monte Carlo (MC) simulation method. Our results show that uncertainties in the studied domain are highly localized. Model simplifications, such as disregarding the groundwater flow, lead to overly confident predictions, that is, predictions that are both less accurate and uncertain than those of the more complex model. We find that for the same number of uncertain parameters, increasing the model resolution reduces uncertainty in the model predictions (and increases the MC method's computational cost). We employ the multilevel MC method to reduce the cost of estimating uncertainty in a high‐resolution ICPR model. Finally, we use the ensemble estimates of the mean and covariance of the flood depth for real‐time flood depth forecasting using the physics‐informed Gaussian process regression method. We show that even with few measurements, the proposed framework results in a more accurate forecast than that provided by the mean prediction of the ICPR model.

A Multilevel Monte Carlo Approach for a Stochastic Optimal Control Problem Based on the Gradient Projection Method

A multilevel Monte Carlo (MLMC) method is applied to simulate a stochastic optimal problem based on the gradient projection method. In the numerical simulation of the stochastic optimal control problem, the approximation of expected value is involved, and the MLMC method is used to address it. The computational cost of the MLMC method and the convergence analysis of the MLMC gradient projection algorithm are presented. Two numerical examples are carried out to verify the effectiveness of our method.

Multi-fidelity microstructure-induced uncertainty quantification by advanced Monte Carlo methods

Quantifying uncertainty associated with the microstructure variation of a material can be a computationally daunting task, especially when dealing with advanced constitutive models and fine mesh resolutions in the crystal plasticity finite element method (CPFEM). Numerous studies have been conducted regarding the sensitivity of material properties and performance to the mesh resolution and choice of constitutive model. However, a unified approach that accounts for various fidelity parameters, such as mesh resolutions, integration time-steps and constitutive models simultaneously is currently lacking. This paper proposes a novel uncertainty quantification (UQ) approach for computing the properties and performance of homogenized materials using CPFEM, that exploits a hierarchy of approximations with different levels of fidelity. In particular, we illustrate how multi-level sampling methods, such as multi-level Monte Carlo (MLMC) and multi-index Monte Carlo (MIMC), can be applied to assess the impact of variations in the microstructure of polycrystalline materials on the predictions of homogenized materials properties. We show that by adaptively exploiting the fidelity hierarchy, we can significantly reduce the number of microstructures required to reach a certain prescribed accuracy. Finally, we show how our approach can be extended to a multi-fidelity framework, where we allow the underlying constitutive model to be chosen from either a phenomenological plasticity model or a dislocation-density-based model.

Quantify uncertainty by estimating the probability density function of the output of interest using MLMC based Bayes method

<p style='text-indent:20px;'>In uncertainty quantification, the quantity of interest is usually the statistics of the space and/or time integration of system solution. In order to reduce the computational cost, a Bayes estimator based on multilevel Monte Carlo (MLMC) is introduced in this paper. The cumulative distribution function of the output of interest, that is, the expectation of the indicator function, is estimated by MLMC method instead of the classic Monte Carlo simulation. Then, combined with the corresponding probability density function, the quantity of interest is obtained by using some specific quadrature rules. In addition, the smoothing of indicator function and Latin hypercube sampling are used to accelerate the reduction of variance. An elliptic stochastic partial differential equation is used to provide a research context for this model. Numerical experiments are performed to verify the advantage of computational reduction and accuracy improvement of our MLMC-Bayes method.</p>

QUANTIFYING UNCERTAIN SYSTEM OUTPUTS VIA THE MULTI-LEVEL MONTE CARLO METHOD-DISTRIBUTION AND ROBUSTNESS MEASURES

In this work, we consider the problem of estimating the probability distribution, the quantile or the conditional expectation above the quantile, the so called conditional-value-at-risk (CVaR), of output quantities of complex random differential models by the Multilevel Monte Carlo (MLMC) method. We follow an approach that recasts the estimation of the above quantities to the computation of suitable parametric expectations. In this work, we present novel computable error estimators for the estimation of such quantities, which are then used to optimally tune the MLMC hierarchy in a continuation type adaptive algorithm. We demonstrate the efficiency and robustness of our adaptive continuation-MLMC in an array of numerical test cases.

Unbiased Optimal Stopping via the MUSE

We propose a new unbiased estimator for estimating the utility of the optimal stopping problem. The MUSE, short for ‘Multilevel Unbiased Stopping Estimator’, constructs the unbiased Multilevel Monte Carlo (MLMC) estimator at every stage of the optimal stopping problem in a backward recursive way. In contrast to traditional sequential methods, the MUSE can be implemented in parallel. We prove the MUSE has finite variance, finite computational complexity, and achieves ɛ-accuracy with O(1/ɛ2) computational cost under mild conditions. We demonstrate MUSE empirically in an option pricing problem involving a high-dimensional input and the use of many parallel processors.