Abstract. In this article, we derive a series expansion of the multivariate normal probability integrals based on Fourier series. The basic idea is to transform the limits of each integral from hi to ∞ to be from −∞ to ∞ by multiplying the integrand by a periodic square wave that approximates the domain of the integral. This square wave is expressed by its Fourier series expansion. Then a Cholesky decomposition of the covariance matrix is applied to transform the integrand to a simple one that can be easily evaluated. The resultant formula has a simple pattern that is expressed as multiple series expansion of trigonometric and exponential functions.