  We apply the holonomic gradient method (HGM) introduced by [9] to the
calculation of orthant probabilities of multivariate normal distribution. The
holonomic gradient method applied to orthant probabilities is found to be a
variant of Plackett's recurrence relation ([14]). However an implementation of
the method yields recurrence relations more suitable for numerical computation
than Plackett's recurrence relation. We derive some theoretical results on the
holonomic system for the orthant probabilities. These results show that
multivariate normal orthant probabilities possess some remarkable properties
from the viewpoint of holonomic systems. Finally we show that numerical
performance of our method is comparable or superior compared to existing
methods.
