High-dimensional Integrals Research Articles

Pricing Path-Dependent Payoff Using Quadrature on Wiener Space

Solution to stochastic differential equation that defines the dynamics of asset prices is interpreted as a function of underlying Brownian motion when permissible, and the expectation values of various functionals of the solution – therefore functionals of Brownian motion (Wiener functionals) – are estimated from numerical approximations to high-dimensional Gaussian integrals using quadrature. The simple method is exemplified for a few standard asset price dynamics used in practice for structured and contingent claim pricing, such as local volatility, Heston stochastic volatility, SABR, and LIBOR market model (LMM/BGM). The quadrature in these concrete implementations is performed in m-dimensional space where m is between 10 and 30 depending on the specific model or a payoff structure. The methodology appears to be a broadly viable alternative to Monte Carlo for short- to medium-dated financial products, with potential uses outside those examined in this paper.

Importance Nested Sampling and the MultiNest Algorithm

Bayesian inference involves two main computational challenges. First, in estimating the parameters of some model for the data, the posterior distribution may well be highly multi-modal: a regime in which the convergence to stationarity of traditional Markov Chain Monte Carlo (MCMC) techniques becomes incredibly slow. Second, in selecting between a set of competing models the necessary estimation of the Bayesian evidence for each is, by definition, a (possibly high-dimensional) integration over the entire parameter space; again this can be a daunting computational task, although new Monte Carlo (MC) integration algorithms offer solutions of ever increasing efficiency. Nested sampling (NS) is one such contemporary MC strategy targeted at calculation of the Bayesian evidence, but which also enables posterior inference as a by-product, thereby allowing simultaneous parameter estimation and model selection. The widely-used MultiNest algorithm presents a particularly efficient implementation of the NS technique for multi-modal posteriors. In this paper we discuss importance nested sampling (INS), an alternative summation of the MultiNest draws, which can calculate the Bayesian evidence at up to an order of magnitude higher accuracy than `vanilla' NS with no change in the way MultiNest explores the parameter space. This is accomplished by treating as a (pseudo-)importance sample the totality of points collected by MultiNest, including those previously discarded under the constrained likelihood sampling of the NS algorithm. We apply this technique to several challenging test problems and compare the accuracy of Bayesian evidences obtained with INS against those from vanilla NS.

Derandomised lattice rules for high dimensional integration

We seek shifted lattice rules that are good for high dimensional integration over the unit cube in the setting of an unanchored weighted Sobolev space of functions with square-integrable mixed first derivatives. Many existing studies rely on random shifting of the lattice, whereas here we work with lattice rules with a deterministic shift. Specifically, we consider 'half-shifted' rules in which each component of the shift is an odd multiple of \(1/(2N)\) where \(N\) is the number of points in the lattice. By applying the principle that there is always at least one choice as good as the average, we show that for a given generating vector there exists a half-shifted rule whose squared worst-case error differs from the shift-averaged squared worst-case error by a term of only order \({1/N^2}\). We carry out numerical experiments where the generating vector is chosen component-by-component (CBC), as for randomly shifted lattices, and where the shift is chosen by a new `CBC for shift' algorithm. The numerical results are encouraging. 
 
 References J. Dick, F. Y. Kuo, and I. H. Sloan. High-dimensional integration: The quasi-Monte Carlo way. Acta Numer., 22:133–288, 2013. doi:10.1017/S0962492913000044. J. Dick, D. Nuyens, and F. Pillichshammer. Lattice rules for nonperiodic smooth integrands. Numer. Math., 126(2):259–291, 2014. doi:10.1007/s00211-013-0566-0. T. Goda, K. Suzuki, and T. Yoshiki. Lattice rules in non-periodic subspaces of sobolev spaces. Numer. Math., 141(2):399–427, 2019. doi:10.1007/s00211-018-1003-1. F. Y. Kuo. Lattice rule generating vectors. URL http://web.maths.unsw.edu.au/ fkuo/lattice/index.html. D. Nuyens and R. Cools. Fast algorithms for component-by-component construction of rank-1 lattice rules in shift-invariant reproducing kernel Hilbert spaces. Math. Comput., 75:903–920, 2006. doi:10.1090/S0025-5718-06-01785-6. I. H. Sloan and S. Joe. Lattice methods for multiple integration. Oxford Science Publications. Clarendon Press and Oxford University Press, 1994. URL https://global.oup.com/academic/product/lattice-methods-for-multiple-integration-9780198534723. I. H. Sloan and H. Wozniakowski. When are quasi-Monte Carlo algorithms efficient for high dimensional integrals? J. Complex., 14(1):1–33, 1998. doi:10.1006/jcom.1997.0463. I. H. Sloan, F. Y. Kuo, and S. Joe. On the step-by-step construction of quasi-Monte Carlo integration rules that achieve strong tractability error bounds in weighted Sobolev spaces. Math. Comput., 71:1609–1641, 2002. doi:10.1090/S0025-5718-02-01420-5.

A Monte Carlo framework for missing wedge restoration and noise removal in cryo-electron tomography

We propose a statistical method to address an important issue in cryo-electron tomography image analysis: reduction of a high amount of noise and artifacts due to the presence of a missing wedge (MW) in the spectral domain. The method takes as an input a 3D tomogram derived from limited-angle tomography, and gives as an output a 3D denoised and artifact compensated volume. The artifact compensation is achieved by filling up the MW with meaningful information. To address this inverse problem, we compute a Minimum Mean Square Error (MMSE) estimator of the uncorrupted image. The underlying high-dimensional integral is computed by applying a dedicated Markov Chain Monte-Carlo (MCMC) sampling procedure based on the Metropolis-Hasting (MH) algorithm. The proposed MWR (Missing Wedge Restoration) algorithm can be used to enhance visualization or as a pre-processing step for image analysis, including segmentation and classification of macromolecules. Results are presented for both synthetic data and real 3D cryo-electron images.

Parallel cross interpolation for high-precision calculation of high-dimensional integrals

We propose a parallel version of the cross interpolation algorithm and apply it to calculate high-dimensional integrals motivated by Ising model in quantum physics. In contrast to mainstream approaches, such as Monte Carlo and quasi Monte Carlo, the samples calculated by our algorithm are neither random nor form a regular lattice. Instead we calculate the given function along individual dimensions (modes) and use these values to reconstruct its behaviour in the whole domain. The positions of the calculated univariate fibres are chosen adaptively for the given function. The required evaluations can be executed in parallel along each mode (variable) and over all modes.To demonstrate the efficiency of the proposed method, we apply it to compute high-dimensional Ising susceptibility integrals, arising from asymptotic expansions for the spontaneous magnetisation in two-dimensional Ising model of ferromagnetism. We observe strong superlinear convergence of the proposed method, while the MC and qMC algorithms converge sublinearly. Using multiple precision arithmetic, we also observe exponential convergence of the proposed algorithm. Combining high-order convergence, almost perfect scalability up to hundreds of processes, and the same flexibility as MC and qMC, the proposed algorithm can be a new method of choice for problems involving high-dimensional integration, e.g. in statistics, probability, and quantum physics.

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.

Restricted maximum likelihood estimation of a censored random effects panel regression model

Panel data sets have been developed in various areas, and many recent studies have analyzed panel, or longitudinal data sets. Maximum likelihood (ML) may be the most common statistical method for analyzing panel data models; however, the inference based on the ML estimate will have an inflated Type I error because the ML method tends to give a downwardly biased estimate of variance components when the sample size is small. The under estimation could be severe when data is incomplete. This paper proposes the restricted maximum likelihood (REML) method for a random effects panel data model with a censored dependent variable. Note that the likelihood function of the model is complex in that it includes a multidimensional integral. Many authors proposed to use integral approximation methods for the computation of likelihood function; however, it is well known that integral approximation methods are inadequate for high dimensional integrals in practice. This paper introduces to use the moments of truncated multivariate normal random vector for the calculation of multidimensional integral. In addition, a proper asymptotic standard error of REML estimate is given.

Improved cross entropy-based importance sampling with a flexible mixture model

The probability of a rare event or failure event is defined through a potentially high-dimensional integral, whose integration domain is often only known point-wise in terms of the outcome of a numerical model. The probability of failure can be estimated efficiently through importance sampling (IS), provided that an effective IS density is chosen. The cross entropy (CE) method is an adaptive sampling approach that determines the IS density through minimizing the Kullback–Leibler divergence between the theoretically optimal IS density and a chosen parametric family of distributions. We propose an improved version of the classical CE method that introduces a smooth transition to make better use of the samples from intermediate sampling levels for fitting the sought IS density. The improved CE method is combined with a novel flexible parametric distribution model that is able to handle low- and high-dimensional problems as well as problems with multimodal failure domains. A set of numerical examples demonstrate that the proposed approach performs consistently better than the classical CE method in various problem settings.

An Automated Prior Robustness Analysis in Bayesian Model Comparison

[enter Abstract BThe marginal likelihood is the gold standard for Bayesian model comparison although it is well-known that the value of marginal likelihood could be sensitive to the choice of prior hyperparameters. Most models require computationally intense simulation-based methods to evaluate the typically high-dimensional integral of the marginal likelihood expression. Hence, despite the recognition that prior sensitivity analysis is important in this context, it is rarely done in practice. In this paper we develop efficient and feasible methods to compute the sensitivities of marginal likelihood, obtained via two common simulation-based methods, with respect to any prior hyperparameter alongside the MCMC estimation algorithm. Our approach builds on Automatic Differentiation (AD), which has only recently been introduced to the more computationally intensive setting of Markov chain Monte Carlo simulation. We illustrate our approach with two empirical applications in the context of widely used multivariate time series models.

Dependence modeling for multi-type recurrent events via copulas.

When several types of recurrent events may arise, interest often lies in marginal modeling and studying the nature of the dependence structure. In this paper, we propose a multivariate mixed-Poisson model with the dependence between events accommodated by type-specific random effects which are associated through use of a Gaussian copula. Such models retain marginal features with a simple interpretation, reflect the heterogeneity in risk for each type of event, and provide insight into the dependence between the different types of events. Semiparametric inference is proposed based on composite likelihood to avoid high dimensional integration. An application to a study of nutritional supplements in malnourished children is given in which the goal is to evaluate the reduction in the rate of several different kinds of infection.

Decrease of the mean of the quasi-random integration error

The article is devoted to the study of the behavior of the quasi-random integration remainder in the calculation of high-dimensional integrals. As noted in the previous work of the authors, the asymptotic behavior of its decrease, determined by the Koksma-Hlawka inequality, can be used only with a very large number of integration nodes N, which cannot be implemented on modern computers. The article introduces the concept of a mean order of decreasing remainder, which makes it possible to judge its properties with the N values available for realization and to compare various pseudo-random sequences. A number of numerical examples are given. In all cases, it turned out that the Sobol’ sequences in the sense of this criterion are somewhat better than the Holton sequences.

IMPRes-Pro: A high dimensional multiomics integration method for in silico hypothesis generation.

Nowadays, large amounts of omics data have been generated and contributed to increasing knowledge about associated biological mechanisms. A new challenge coming along is how to identify the active pathways and extract useful insights from these data with huge background information and noise. Although biologically meaningful modules can often be detected by many existing informatics tools, it is still hard to interpret or make use of the results towards in silico hypothesis generation and testing. To address this gap, we previously developed the IMPRes (Integrative MultiOmics Pathway Resolution) v 1.0 algorithm, a new step-wise active pathway detection method using a dynamic programming approach. This approach enables the network detection one step at a time, making it easy for researchers to trace the pathways, and leading to more accurate drug design and more effective treatment strategies. In this paper, we present IMPRes-Pro, an enhancement to IMPRes v1.0 by integrating proteomics data along with transcriptomics data and constructing a heterogeneous background network. The evaluation experiment conducted on human primary breast cancer dataset has shown the advantage over the original IMPRes v1.0 method. Furthermore, a case study on human metastatic breast cancer dataset was performed and we have provided several insights regarding the selection of optimal therapy strategy. IMPRes-Pro algorithm and visualization tool is available as a web service at http://digbio.missouri.edu/impres.

Analysis of quasi-Monte Carlo methods for elliptic eigenvalue problems with stochastic coefficients

We consider the forward problem of uncertainty quantification for the generalised Dirichlet eigenvalue problem for a coercive second order partial differential operator with random coefficients, motivated by problems in structural mechanics, photonic crystals and neutron diffusion. The PDE coefficients are assumed to be uniformly bounded random fields, represented as infinite series parametrised by uniformly distributed i.i.d. random variables. The expectation of the fundamental eigenvalue of this problem is computed by (a) truncating the infinite series which define the coefficients; (b) approximating the resulting truncated problem using lowest order conforming finite elements and a sparse matrix eigenvalue solver; and (c) approximating the resulting finite (but high dimensional) integral by a randomly shifted quasi-Monte Carlo lattice rule, with specially chosen generating vector. We prove error estimates for the combined error, which depend on the truncation dimension s, the finite element mesh diameter h, and the number of quasi-Monte Carlo samples N. Under suitable regularity assumptions, our bounds are of the particular form $${\mathcal {O}}(h^2 + N^{-1 + \delta })$$ , where $$\delta > 0$$ is arbitrary and the hidden constant is independent of the truncation dimension, which needs to grow as $$h\rightarrow 0$$ and $$N \rightarrow \infty $$ . As for the analogous PDE source problem, the conditions under which our error bounds hold depend on a parameter $$p \in (0, 1)$$ representing the summability of the terms in the series expansions of the coefficients. Although the eigenvalue problem is nonlinear, which means it is generally considered harder than the source problem, in almost all cases ( $$p \ne 1$$ ) we obtain error bounds that converge at the same rate as the corresponding rate for the source problem. The proof involves a detailed study of the regularity of the fundamental eigenvalue as a function of the random parameters. As a key intermediate result in the analysis, we prove that the spectral gap (between the fundamental and the second eigenvalues) is uniformly positive over all realisations of the random problem.

Maximum likelihood estimation of Gaussian copula models for geostatistical count data

This work investigates the computation of maximum likelihood estimators in Gaussian copula models for geostatistical count data. This is a computationally challenging task because the likelihood function is only expressible as a high dimensional multivariate normal integral. Two previously proposed Monte Carlo methods are reviewed, the Genz–Bretz and Geweke–Hajivassiliou–Keane simulators, and a new method is investigated. The new method is based on the so–called data cloning algorithm, which uses Markov chain Monte Carlo algorithms to approximate maximum likelihood estimators and their (asymptotic) variances in models with computationally challenging likelihoods. A simulation study is carried out to compare the statistical and computational efficiencies of the three methods. It is found that the three methods have similar statistical properties, but the Geweke–Hajivassiliou–Keane simulator requires the least computational effort. Hence, this is the method we recommend. A data analysis of Lansing Woods tree counts is used to illustrate the methods.

Composite Likelihood Estimation of an Autoregressive Panel Probit Model with Random Effects

Modeling and estimating autocorrelated discrete data can be challenging. In this paper, we use an autoregressive panel probit model where the autocorrelation in the discrete variable is driven by the autocorrelation in the latent variable. In such a non-linear model, the autocorrelation in an unobserved variable results in an intractable likelihood containing high-dimensional integrals. To tackle this problem, we use composite likelihoods that involve much lower order of integration. However, parameter identification becomes problematic since the information employed in lower dimensional distributions may not be rich enough for identification. Therefore, we characterize types of composite likelihoods that are valid for this model and study conditions under which the parameters can be identified. Moreover, we provide consistency and asymptotic normality results of the pairwise composite likelihood estimator and conduct Monte Carlo studies to assess its finite-sample performances. Finally, we apply our method to analyze credit ratings. The results indicate a significant improvement in the estimated probabilities for rating transitions compared with static models.

A continuous analogue of the tensor-train decomposition

We develop new approximation algorithms and data structures for representing and computing with multivariate functions using the functional tensor-train (FT), a continuous extension of the tensor-train (TT) decomposition. The FT represents functions using a tensor-train ansatz by replacing the three-dimensional TT cores with univariate matrix-valued functions. The main contribution of this paper is a framework to compute the FT that employs adaptive approximations of univariate fibers, and that is not tied to any tensorized discretization. The algorithm can be coupled with any univariate linear or nonlinear approximation procedure. We demonstrate that this approach can generate multivariate function approximations that are several orders of magnitude more accurate, for the same cost, than those based on the conventional approach of compressing the coefficient tensor of a tensor-product basis. Our approach is in the spirit of other continuous computation packages such as Chebfun, and yields an algorithm which requires the computation of “continuous” matrix factorizations such as the LU and QR decompositions of vector-valued functions. To support these developments, we describe continuous versions of an approximate maximum-volume cross approximation algorithm and of a rounding algorithm that re-approximates an FT by one of lower ranks. We demonstrate that our technique improves accuracy and robustness, compared to TT and quantics-TT approaches with fixed parameterizations, of high-dimensional integration, differentiation, and approximation of functions with local features such as discontinuities and other nonlinearities.

Stable and efficient cubature rules by metaheuristic optimization with application to Kalman filtering

The numerical approximation of Gaussian integrals is of great relevance for a wide variety of scientific applications. In recent decades, numerous efficient cubature rules of high polynomial exactness for the case of high-dimensional Gaussian integrals have been proposed. In addition to the question of efficiency, however, the question of stability, which depends on the influence of negative weights, is also of crucial importance. For certain degrees of polynomial exactness no stable rules are known. Therefore, a new approach based on metaheuristic optimization is presented which allows the computation of cubature rules which exceed comparable rules known from the literature with respect to stability and efficiency. Especially with regard to iterative numerical calculations, the use of unstable cubature rules over time can lead to a distortion and thus a significant deterioration of the results. To illustrate this fact, the cubature-based Kalman filter is utilized. Based on the results of two simulation studies, the necessity for the use of stable cubature rules is motivated.

A Simple Method for Comparing Complex Models: Bayesian Model Comparison for Hierarchical Multinomial Processing Tree Models Using Warp-III Bridge Sampling

Multinomial processing trees (MPTs) are a popular class of cognitive models for categorical data. Typically, researchers compare several MPTs, each equipped with many parameters, especially when the models are implemented in a hierarchical framework. A Bayesian solution is to compute posterior model probabilities and Bayes factors. Both quantities, however, rely on the marginal likelihood, a high-dimensional integral that cannot be evaluated analytically. In this case study, we show how Warp-III bridge sampling can be used to compute the marginal likelihood for hierarchical MPTs. We illustrate the procedure with two published data sets and demonstrate how Warp-III facilitates Bayesian model averaging.

Coulomb-Volkov S-matrix theory of ionisation

In this paper we outline the recently derived Coulomb-Volkov (CV) S-matrix theory and present the results of calculations of the energy and the momentum spectra of the electrons emitted from ionisation of Ar atoms by an intense long-wavelength laser field. The results are compared with the corresponding calculations using the usual (plane-wave) SFA S-matrix theory, up to the second order of the series. We introduce a modified stationary phase method with a parametric integration to evaluate the high-dimensional integrals that occur in the second order of the ionisation amplitudes. The final-state Coulomb effect is found to be dominant at the lower energies, where the two theories differ greatly. In contrast, at high energies, the two theories produce qualitatively similar (in the logarithmic scale) distributions. It is found that despite the similarity of the positions of the ‘forward rescattering’ (or ‘soft collision’) peaks shown by the second order CV- and the SFA-amplitude, the line-shapes produced by them are qualitatively different. For example, the CV-spectrum shows a peak at 0.1 Up with a characteristic line-shape that is convex on one side and steeply rising on the other side; this is qualitatively similar to the line-shape of the peak observed experimentally earlier at the same energy. In contrast, the SFA-spectrum shows at the same energy a pointed peak with a line-shape that is concave on both sides. The characteristic line-shape of the low energy LES-peak(s) thus appears to suggest itself as a signature of the final-state Coulomb interaction in strong-field ionisation processes.

Deep learning algorithm for data-driven simulation of noisy dynamical system

We present a deep learning model, DE-LSTM, for the simulation of a stochastic process with an underlying nonlinear dynamics. The deep learning model aims to approximate the probability density function of a stochastic process via numerical discretization and the underlying nonlinear dynamics is modeled by the Long Short-Term Memory (LSTM) network. It is shown that, when the numerical discretization is used, the function estimation problem can be solved by a multi-label classification problem. A penalized maximum log likelihood method is proposed to impose a smoothness condition in the prediction of the probability distribution. We show that the time evolution of the probability distribution can be computed by a high-dimensional integration of the transition probability of the LSTM internal states. A Monte Carlo algorithm to approximate the high-dimensional integration is outlined. The behavior of DE-LSTM is thoroughly investigated by using the Ornstein–Uhlenbeck process and noisy observations of nonlinear dynamical systems; Mackey–Glass time series and forced Van der Pol oscillator. It is shown that DE-LSTM makes a good prediction of the probability distribution without assuming any distributional properties of the stochastic process. For a multiple-step forecast of the Mackey–Glass time series, the prediction uncertainty, denoted by the 95% confidence interval, first grows, then dynamically adjusts following the evolution of the system, while in the simulation of the forced Van der Pol oscillator, the prediction uncertainty does not grow in time even for a 3,000-step forecast.