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

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 100× over packages such as and 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.