Doubly blocked variant of Genz’s method for multivariate normal integration

Abstract

High-dimensional multivariate normal (MVN) integration is a computational bottleneck in many statistical applications, particularly in finance and econometrics. This paper presents a modern, doubly blocked variant of Genz’s method for integrating MVN distributions, implemented in the Julia programming language. By leveraging multithreading, cache optimization, and SIMD/BLAS-friendly design, we address the memory bandwidth limitations of traditional implementations. We benchmark our approach against existing solvers in Julia, Fortran, and C++, demonstrating speedups that are often close to or above \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$100\times $$\end{document} over packages such as mvtnorm and tlrmvnmvt for large dimensions in the tested settings. We also compare our method with tiled low-rank approximations and apply it to a high-dimensional Gaussian copula model using real-world financial data. The results show that our implementation enables very fast probability calculations for high-dimensional portfolios, a task which was previously hard to achieve with standard tools.