Multivariate phase relationships are important to characterize and understand numerous physical, biological, and chemical systems, from electromagnetic waves to neural oscillations. These systems exhibit complex spatiotemporal dynamics and intricate interdependencies among their constituent elements. While classical models of multivariate phase relationships, such as the wave equation and Kuramoto model, give theoretical models to describe phenomena, the development of statistical tools for hypothesis testing and inference for multivariate phase relationships in complex systems remains limited. This paper introduces a novel probabilistic modeling framework to characterize multivariate phase relationships, with wave-like phenomena serving as a key example. This approach describes spatial patterns and interactions between oscillators through a pairwise exponential family distribution. Building upon the literature of graphical model inference, including methods like Ising models, graphical lasso, and interaction screening, this work bridges the gap between classical wave dynamics and modern statistical approaches. Efficient inference methods are introduced, leveraging the Chow-Liu algorithm for directed tree approximations and interaction screening for general graphical models. Simulated experiments demonstrate the utility of these methods for uncovering wave properties and sparse interaction structures, highlighting their applicability to diverse scientific domains. This framework establishes a new paradigm for statistical modeling of multivariate phase relationships, providing a powerful toolset for exploring the complexity of these systems.