Abstract Developments in the stability of modern spectrographs have led to extremely precise instrumental radial velocity (RV) measurements. For most stars, the detection limit of planetary companions with these instruments is expected to be dominated by astrophysical noise sources such as starspots. Correlated signals caused by rotationally modulated starspots can obscure or mimic the Doppler shifts induced by even the closest, most massive planets. This is especially true for young, magnetically active stars where stellar activity can cause fluctuation amplitudes of ≳0.1 mag in brightness and ≳100 m s⁻¹ in RV semiamplitudes. Techniques that can mitigate these effects and increase our sensitivity to young planets are critical to improving our understanding of the evolution of planetary systems. Gaussian processes (GPs) have been successfully employed to model and constrain activity signals in individual cases. However, a principled approach of this technique, specifically for the joint modeling of photometry and RVs, has not yet been developed. In this work, we present a GP framework to simultaneously model stellar activity signals in photometry and RVs that can be used to investigate the relationship between both time series. Our method, inspired by the FF ′ framework of Aigrain et al., models spot-driven activity signals as the linear combinations of two independent latent GPs and their time derivatives. We also simulate time series affected by starspots by extending the starry software to incorporate time evolution of stellar features. Using these synthetic data sets, we show that our method can predict spot-driven RV variations with greater accuracy than other GP approaches.