Active galactic nuclei (AGN) are believed to be powered by the accretion of matter around supermassive black holes at the centers of galaxies. The variability of an AGN's brightness over time can reveal important information about the physical properties of the underlying black hole. The temporal variability is believed to follow a stochastic process, often represented as a damped random walk described by a stochastic differential equation (SDE). With upcoming wide-field surveys set to observe 100 million AGN in multiple bandpass filters, there is a need for efficient and automated modeling techniques that can handle the large volume of data. Latent SDEs are well-suited for modeling AGN time series data, as they can explicitly capture the underlying stochastic dynamics. In this work, we modify latent SDEs to jointly reconstruct the unobserved portions of multivariate AGN light curves and infer their physical properties such as the black hole mass. Our model is trained on a realistic physics-based simulation of ten-year AGN light curves, and we demonstrate its ability to fit AGN light curves even in the presence of long seasonal gaps and irregular sampling across different bands, outperforming a multi-output Gaussian process regression baseline.