quasi-Monte Carlo Research Articles

Stability of lattice rules and polynomial lattice rules constructed by the component-by-component algorithm

We study quasi-Monte Carlo (QMC) methods for numerical integration of multivariate functions defined over the high-dimensional unit cube. Lattice rules and polynomial lattice rules, which are special classes of QMC methods, have been intensively studied and the so-called component-by-component (CBC) algorithm has been well-established to construct rules which achieve the almost optimal rate of convergence with good tractability properties for given smoothness and set of weights. Since the CBC algorithm constructs rules for given smoothness and weights, not much is known when such rules are used for function classes with different smoothness and/or weights.In this paper we prove that a lattice rule constructed by the CBC algorithm for the weighted Korobov space with given smoothness and weights achieves the almost optimal rate of convergence with good tractability properties for general classes of smoothness and weights which satisfy some summability conditions. Such a stability result also can be shown for polynomial lattice rules in weighted Walsh spaces. We further give bounds on the weighted star discrepancy and discuss the tractability properties for these QMC rules. The results are comparable to those obtained for Halton, Sobol and Niederreiter sequences.

Faster Monte Carlo estimation of joint models for time-to-event and multivariate longitudinal data

Quasi-Monte Carlo (QMC) methods using quasi-random sequences, as opposed to pseudo-random samples, are proposed for use in the joint modelling of time-to-event and multivariate longitudinal data. The QMC integration framework extends the Monte Carlo Expectation Maximisation approaches that are commonly adopted, namely using ordinary and antithetic variates. The motivation of QMC integration is to increase the convergence speed by using nodes that are scattered more uniformly. Through simulation, estimates and computational times are compared and this is followed with an application to a clinical dataset. There is a distinct speed advantage in using QMC methods for small sample sizes and QMC is comparable to the antithetic MC method for moderate sample sizes. The new method is available in an updated version of the R package joineRML.

Primal–dual quasi-Monte Carlo simulation with dimension reduction for pricing American options

The pricing of American options is one of the most challenging problems in financial engineering due to the involved optimal stopping time problem, which can be solved by using dynamic programming (DP). But applying DP is not always practical, especially when the state space is high dimensional. However, the curse of dimensionality can be overcome by Monte Carlo (MC) simulation. We can get lower and upper bounds by MC to ensure that the true price falls into a valid confidence interval. During the recent decades, progress has been made in using MC simulation to obtain both the lower bound by least-squares Monte Carlo method (LSM) and the upper bound by duality approach. However, there are few works on pricing American options using quasi-Monte Carlo (QMC) methods, especially to compute the upper bound. For comparing the sample variances and standard errors in the numerical experiments, randomized QMC (RQMC) methods are usually used. In this paper, we propose to use RQMC to replace MC simulation to compute both the lower bound (by the LSM) and the upper bound (by the duality approach). Moreover, we propose to use dimension reduction techniques, such as the Brownian bridge, principal component analysis, linear transformation and the gradients based principle component analysis. We perform numerical experiments on American–Asian options and American max-call options under the Black–Scholes model and the variance gamma model, in which the options have the path-dependent feature or are written on multiple underlying assets. We find that RQMC in combination with dimension reduction techniques can significantly increase the efficiency in computing both the lower and upper bounds, resulting in better estimates and tighter confidence intervals of the true price than pure MC simulation.

P-Refined Multilevel Quasi-Monte Carlo for Galerkin Finite Element Methods with Applications in Civil Engineering

Civil engineering applications are often characterized by a large uncertainty on the material parameters. Discretization of the underlying equations is typically done by means of the Galerkin Finite Element method. The uncertain material parameter can be expressed as a random field represented by, for example, a Karhunen–Loève expansion. Computation of the stochastic responses, i.e., the expected value and variance of a chosen quantity of interest, remains very costly, even when state-of-the-art Multilevel Monte Carlo (MLMC) is used. A significant cost reduction can be achieved by using a recently developed multilevel method: p-refined Multilevel Quasi-Monte Carlo (p-MLQMC). This method is based on the idea of variance reduction by employing a hierarchical discretization of the problem based on a p-refinement scheme. It is combined with a rank-1 Quasi-Monte Carlo (QMC) lattice rule, which yields faster convergence compared to the use of random Monte Carlo points. In this work, we developed algorithms for the p-MLQMC method for two dimensional problems. The p-MLQMC method is first benchmarked on an academic beam problem. Finally, we use our algorithm for the assessment of the stability of slopes, a problem that arises in geotechnical engineering, and typically suffers from large parameter uncertainty. For both considered problems, we observe a very significant reduction in the amount of computational work with respect to MLMC.

Hierarchical adaptive sparse grids and quasi-Monte Carlo for option pricing under the rough Bergomi model

The rough Bergomi (rBergomi) model, introduced recently in Bayer et al. [Pricing under rough volatility. Quant. Finance, 2016, 16(6), 887–904], is a promising rough volatility model in quantitative finance. It is a parsimonious model depending on only three parameters, and yet remarkably fits empirical implied volatility surfaces. In the absence of analytical European option pricing methods for the model, and due to the non-Markovian nature of the fractional driver, the prevalent option is to use the Monte Carlo (MC) simulation for pricing. Despite recent advances in the MC method in this context, pricing under the rBergomi model is still a time-consuming task. To overcome this issue, we have designed a novel, hierarchical approach, based on: (i) adaptive sparse grids quadrature (ASGQ), and (ii) quasi-Monte Carlo (QMC). Both techniques are coupled with a Brownian bridge construction and a Richardson extrapolation on the weak error. By uncovering the available regularity, our hierarchical methods demonstrate substantial computational gains with respect to the standard MC method. They reach a sufficiently small relative error tolerance in the price estimates across different parameter constellations, even for very small values of the Hurst parameter. Our work opens a new research direction in this field, i.e. to investigate the performance of methods other than Monte Carlo for pricing and calibrating under the rBergomi model.

An Integrated Quasi-Monte Carlo Method for Handling High Dimensional Problems with Discontinuities in Financial Engineering

Quasi-Monte Carlo (QMC) method is a useful numerical tool for pricing and hedging of complex financial derivatives. These problems are usually of high dimensionality and discontinuities. The two factors may significantly deteriorate the performance of the QMC method. This paper develops an integrated method that overcomes the challenges of the high dimensionality and discontinuities concurrently. For this purpose, a smoothing method is proposed to remove the discontinuities for some typical functions arising from financial engineering. To make the smoothing method applicable for more general functions, a new path generation method is designed for simulating the paths of the underlying assets such that the resulting function has the required form. The new path generation method has an additional power to reduce the effective dimension of the target function. Our proposed method caters for a large variety of model specifications, including the Black-Scholes, exponential normal inverse Gaussian L\'evy, and Heston models. Numerical experiments dealing with these models show that in the QMC setting the proposed smoothing method in combination with the new path generation method can lead to a dramatic variance reduction for pricing exotic options with discontinuous payoffs and for calculating options' Greeks. The investigation on the effective dimension and the related characteristics explains the significant enhancement of the combined procedure.

Deep learning observables in computational fluid dynamics

Many large scale problems in computational fluid dynamics such as uncertainty quantification, Bayesian inversion, data assimilation and PDE constrained optimization are considered very challenging computationally as they require a large number of expensive (forward) numerical solutions of the corresponding PDEs. We propose a machine learning algorithm, based on deep artificial neural networks, that predicts the underlying input parameters to observable map from a few training samples (computed realizations of this map). By a judicious combination of theoretical arguments and empirical observations, we find suitable network architectures and training hyperparameters that result in robust and efficient neural network approximations of the parameters to observable map. Numerical experiments are presented to demonstrate low prediction errors for the trained network networks, even when the network has been trained with a few samples, at a computational cost which is several orders of magnitude lower than the underlying PDE solver.Moreover, we combine the proposed deep learning algorithm with Monte Carlo (MC) and Quasi-Monte Carlo (QMC) methods to efficiently compute uncertainty propagation for nonlinear PDEs. Under the assumption that the underlying neural networks generalize well, we prove that the deep learning MC and QMC algorithms are guaranteed to be faster than the baseline (quasi-) Monte Carlo methods. Numerical experiments demonstrating one to two orders of magnitude speed up over baseline QMC and MC algorithms, for the intricate problem of computing probability distributions of the observable, are also presented.

Stochastic economic dispatch considering the dependence of multiple wind farms using multivariate Gaussian kernel copula

Wind farms usually gather in areas with abundant wind resources. Therefore, spatial dependence of wind speeds among nearby wind farms should be taken into account when modeling a power system with large-scale wind power penetration. This paper proposes a novel non-parametric copula method, multivariate Gaussian kernel copula (MGKC), to describe the dependence structure of wind speeds among multiple wind farms. Wind speed scenarios considering the dependence among different wind farms are sampled from the MGKC by quasi-Monte Carlo (QMC) method, so as to solve the stochastic economic dispatch (SED) problem, for which an improved mean-variance (MV) model is established, which targets at minimizing the expectation and risk of fuel cost simultaneously. In this model, confidence interval is applied in the wind speed to obtain more practical dispatch solutions by excluding extreme scenarios, for which the quantile-copula is proposed to construct the confidence interval constraint. Simulation studies are carried out on a modified IEEE 30-bus power system with wind farms integrated in two areas, and the results prove the superiority of the MGKC in formulating the dependence among different wind farms and the superiority of the improved MV model based on quantile-copula in determining a better dispatch solution.

A quasi-Monte Carlo solver for thermal radiation in participating media

The Monte Carlo (MC) method is the most accurate method for resolving radiative heat transfer in participating media. However, it is also computationally prohibitive in large-scale simulations. To alleviate this, this study proposes a quasi-Monte Carlo (QMC) method for thermal radiation in participating media with a focus on combustion-related problems. The QMC method employs low-discrepancy sequences (LDS) in place of the traditional random numbers. Three different low-discrepancy sequences – Sobol, Halton, and Niederreiter – were examined as part of this work. The developed QMC method was first validated against analytical solutions of radiative heat transfer in several one-dimensional configurations. Then it was extended to three-dimensional practical combustion configurations. The results from QMC and traditional Monte Carlo are compared against benchmark solutions for each cases. It is shown that the error of the predicted radiation field from QMC is lower than an equivalent MC simulation. The computational cost of QMC was also found lower than MC due to the avoidance of requirement of several statistical runs for traditional Monte Carlo methods alongside achieving the reduction in error. In conclusion, significant improvements in computational costs and accuracy seen in the QMC method makes it an attractive alternative to traditional Monte Carlo methods in high-fidelity simulations.

Quasi-Monte Carlo tractability of integration problem in function spaces defined over products of balls

Recently, quasi-Monte Carlo (QMC) rules for multivariate integration in some weighted Sobolev spaces of functions defined over unit cubes [Formula: see text], products of [Formula: see text] copies of the simplex [Formula: see text] and products of [Formula: see text] copies of the unit sphere [Formula: see text] have been well-studied in the literature. In this paper, we consider QMC tractability of integrals of functions defined over product of [Formula: see text] copies of the ball [Formula: see text]. The space is a tensor product of [Formula: see text] reproducing kernel Hilbert spaces defined by positive and uniformly bounded weights [Formula: see text] for [Formula: see text]. We obtain matching necessary and sufficient conditions in terms of this weights for various notions of QMC tractability, including strong polynomial tractability, polynomial tractability, quasi-polynomial tractability, uniformly weak tractability and [Formula: see text]-weak tractability.

Quasi-Monte Carlo technique in global sensitivity analysis of wind resource assessment with a study on UAE

The present paper bridges mathematical modeling and wind resource assessment (WRA). Sensitivity analysis (SA) links portions of output variance to the variance in each input variable. Global SA (GSA) explores inputs globally. One-at-a-time SA is dominating in WRA, while GSA is often overlooked. Compared to traditional methods, GSA offers potential improvement by the means of the quasi-Monte Carlo (QMC) technique with its elaborate sampling designs enabling faster convergence. The main novelty of this work is twofold: the use of QMC in WRA and Sobol method as a variance-based GSA method in WRA. This paper is among a few using GSA in WRA. Two case studies were conducted. One shows that QMC with sampling based on Sobol sampling outperforms Latin hypercube sampling and pseudorandom sampling. It also shows that in terms of accuracy of results, the brute-force calculation of Sobol sensitivity indices (SI) should be used whenever the model allows it; otherwise, SI can be estimated. Another case study demonstrates a valid GSA study for WRA at Masdar City, United Arab Emirates. The results suggest that the influence of the variance in Weibull parameters on annual energy production (AEP) might be overestimated, while found to be responsible for 2% of AEP uncertainty, and the influence of the variance in air density might be overlooked, while found to account for 94%. WRA studies would benefit greatly from using the QMC and Sobol method. The Sobol method is a universal GSA method, providing valid results for nonlinear models typical for WRA, and QMC provides global scalability, model independence, and flexibility in uncertainty quantification.

Quasi-Monte Carlo method for solving Fredholm equations

Abstract A Monte Carlo method used for the estimation of convergent von Neumann series solutions of a Fredholm equation of second kind is considered. The sum z ( d ) ⁢ ( x ) {z^{(d)}(x)} of d initial terms of the von Neumann series estimating the solution z ⁢ ( x ) {z(x)} of the equation is represented as a d-dimensional integral over the unit cube H d {H_{d}} . This note presents three examples calculating z ( d ) ⁢ ( x ) {z^{(d)}(x)} for different kernels with norms ∥ K ∥ < 1 {\lVert K\rVert<1} . We found that z ( d ) ⁢ ( x ) {z^{(d)}(x)} calculated using a quasi-Monte Carlo (QMC) method converges significantly faster than the corresponding Monte Carlo (MC) estimates in the entire range of ∥ K ∥ {\lVert K\rVert} values. We also found that the average dimension d ^ {\hat{d}} of the integrand in all our examples is small, less than 3. We suggest that the average dimensions d ^ {\hat{d}} of our d-dimensional integrands are bounded as d → ∞ {d\to\infty} .

A computational investigation of the optimal Halton sequence in QMC applications

Abstract We propose the use of randomized (scrambled) quasirandom sequences for the purpose of providing practical error estimates for quasi-Monte Carlo (QMC) applications. One popular quasirandom sequence among practitioners is the Halton sequence. However, Halton subsequences have correlation problems in their highest dimensions, and so using this sequence for high-dimensional integrals dramatically affects the accuracy of QMC. Consequently, QMC studies have previously proposed several scrambling methods; however, to varying degrees, scrambled versions of Halton sequences still suffer from the correlation problem as manifested in two-dimensional projections. This paper proposes a modified Halton sequence (MHalton), created using a linear digital scrambling method, which finds the optimal multiplier for the Halton sequence in the linear scrambling space. In order to generate better uniformity of distributed sequences, we have chosen strong MHalton multipliers up to 360 dimensions. The proposed multipliers have been tested and proved to be stronger than several sets of multipliers used in other known scrambling methods. To compare the quality of our proposed scrambled MHalton sequences with others, we have performed several extensive computational tests that use L 2 {L_{2}} -discrepancy and high-dimensional integration tests. Moreover, we have tested MHalton sequences on Mortgage-backed security (MBS), which is one of the most widely used applications in finance. We have tested our proposed MHalton sequence numerically and empirically, and they show optimal results in QMC applications. These confirm the efficiency and safety of our proposed MHalton over scrambling sequences previously used in QMC applications.

Parametric sensitivity analysis of SEL-421 distance relay algorithms used in compensated line

This paper presents the use of Global Sensitivity Analysis (GSA) methods to analyze the sensitivity level of the performance of the impedance measurement algorithm of SEL-421 distance relay protection when used to perform transmission line protection with compensation circuit. This GSA technique was developed to see the dominance effect of system uncertainty (or factors). Two techniques that might be used for sensitivity analysis using GSA are Quasi-Monte Carlo (QMC) and data sampling with Sobol technique. Experimental is done with DIgSILENT software to do modeling of the power system, simulation, and calculation of impedance and algorithm error. The testing is carried out automatically with the help of the developed algorithm using DPL (DIgSILENT Programming Language). SIMLAB software, in this case, (with Sobol technique) is applied to generate some non-linear factor data and analyze the output of variations in impedance measurement errors.

Quasi-Monte Carlo based uncertainty analysis: Sampling efficiency and error estimation in engineering applications

In this paper, the potential benefits of quasi-Monte Carlo (QMC) methods for uncertainty propagation are assessed via two applications: a numerical case study and a realistic and complex building physical case study. The sampling efficiency of four quasi-Monte Carlo sampling strategies — Optimized Latin hypercube, Sobol’ sequence, Niederreiter–Xing sequence and lattice sequence — are quantified and the errors of these quasi-Monte Carlo methods are estimated. In addition, for getting a better understanding of the potential factors that may influence the performance of quasi-Monte Carlo methods, the effect of the different parameters and the smoothness of the target function for the building physical case study are investigated for two quantities of interest. A sensitivity analysis is performed in which first and total order Sobol’ indices are calculated and a kernel smoother is used to show the effect on the quantity of interest for certain input parameters. The outcomes show that quasi-Monte Carlo methods perform at least as well as the Monte Carlo method and are capable of outperforming the standard Monte Carlo method if the target function is sufficiently smooth and mainly depends on a limited number of parameters.

Accurate and Efficient Explicit Approximations of the Colebrook Flow Friction Equation Based on the Wright ω-Function: Reply to Discussion

This reply gives two corrections of typographical errors in respect to the commented article, and then provides few comments in respect to the discussion and one improved version of the approximation of the Colebrook equation for flow friction, based on the Wright ω-function. Finally, this reply gives an exact explicit version of the Colebrook equation expressed through the Wright ω-function, which does not introduce any additional errors in respect to the original equation. All mentioned approximations are computationally efficient and also very accurate. Results are verified using more than 2 million of Quasi Monte-Carlo samples.

Deterministic Sampling of Expensive Posteriors Using Minimum Energy Designs

Markov chain Monte Carlo (MCMC) methods require a large number of samples to approximate a posterior distribution, which can be costly when the likelihood or prior is expensive to evaluate. The number of samples can be reduced if we can avoid repeated samples and those that are close to each other. This is the idea behind deterministic sampling methods such as quasi-Monte Carlo (QMC). However, the existing QMC methods aim at sampling from a uniform hypercube, which can miss the high probability regions of the posterior distribution and thus the approximation can be poor. Minimum energy design (MinED) is a recently proposed deterministic sampling method, which makes use of the posterior evaluations to obtain a weighted space-filling design in the region of interest. However, the existing implementation of MinED is inefficient because it requires several global optimizations and thus numerous evaluations of the posterior. In this article, we develop an efficient algorithm that can generate MinED samples with few posterior evaluations. We also make several improvements to the MinED criterion to make it perform better in high dimensions. The advantages of MinED over MCMC and QMC are illustrated using an example of calibrating a friction drilling process.

A gPC-intrusive Monte-Carlo scheme for the resolution of the uncertain linear Boltzmann equation

In this paper, we are interested in the resolution of the time-dependent problem of particle transport in a media whose composition is uncertain. The most common resolution strategy consists of running, at prescribed points in uncertain space (at experimental designs points), a simulation device as a black-box. This ensures performing the uncertainty propagation resolution (i.e. resolution of SPDE). This kind of strategy is commonly called non-intrusive. The non-intrusive resolution can be carried with Monte-Carlo, quasi Monte-Carlo, generalized Polynomial Chaos etc. The latter is of interest in this paper for its fast convergence rate. First, we go over and illustrate the main drawbacks of the non-intrusive (gPC or otherwise) uncertainty propagation resolution for the linear Boltzmann equation in a simplified configuration. We then build a new gPC based MC scheme solving intrusively the uncertain counterpart of the problem. The paper ends with some numerical examples.

On the prediction of pseudo relative permeability curves: meta-heuristics versus Quasi-Monte Carlo

This article reports the first application of the Quasi-Monte Carlo (QMC) method for estimation of the pseudo relative permeability curves. In this regards, the performance of several meta-heuristics algorithms have also been compared versus QMC, including the Genetic Algorithm (GA), Particle Swarm Optimization (PSO), and the Artificial Bee Colony (ABC). The mechanism of minimizing the objective-function has been studied, for each method. The QMC has outperformed its counterparts in terms of accuracy and efficiently sweeping the entire search domain. Nevertheless, its computational time requirement is obtained in excess to the meta-heuristics algorithms.

Richardson Extrapolation of Polynomial Lattice Rules

We study multivariate numerical integration of smooth functions in weighted Sobolev spaces with dominating mixed smoothness $\alpha\geq 2$ defined over the $s$-dimensional unit cube. We propose a new quasi-Monte Carlo (QMC)-based quadrature rule, named \emph{extrapolated polynomial lattice rule}, which achieves the almost optimal rate of convergence. Extrapolated polynomial lattice rules are constructed in two steps: i) construction of classical polynomial lattice rules over $\mathbb{F}_b$ with $\alpha$ consecutive sizes of nodes, $b^{m-\alpha+1},\ldots,b^{m}$, and ii) recursive application of Richardson extrapolation to a chain of $\alpha$ approximate values of the integral obtained by consecutive polynomial lattice rules. We prove the existence of good extrapolated polynomial lattice rules achieving the almost optimal order of convergence of the worst-case error in Sobolev spaces with general weights. Then, by restricting to product weights, we show that such good extrapolated polynomial lattice rules can be constructed by the fast component-by-component algorithm under a computable quality criterion. The required total construction cost is of order $(s+\alpha)N\log N$, which improves the currently known result for interlaced polynomial lattice rule, that is of order $s\alpha N\log N$. We also study the dependence of the worst-case error bound on the dimension. A big advantage of our method compared to interlaced polynomial lattice rules is that the fast QMC matrix vector method can be used in this setting, while still achieving the same rate of convergence. Such a method was previously not known. Numerical experiments for test integrands support our theoretical result.