The inversion of geodetic data to estimate the position and depth of an inflating or contracting magmatic body is inherently non-unique. To help constrain a model's geometry, a priori information is often used to fix the depth and position of the source prior to inversion. This is commonly done based on independent information from other datasets such as Magnetotellurics (MT) and seismology which are able to provide an image of the subsurface structure. Here we use geodetic (InSAR and GNSS), seismological and MT data acquired over the Okataina Caldera, New Zealand, and perform a joint inversion where MT are directly included in the inversion of surface deformation data. InSAR and GNSS data acquired over the last ~20 years have shown subsidence across the Okataina caldera and surrounding areas of up to ~15 mm/yr. Assuming the subsidence is largely being driven by the ongoing cooling and contraction of melt, we estimate the distribution of contraction beneath the region. With the inclusion of MT data in the inversion, we discover two important features 1) contraction is largely focussed beneath the two geologically active linear vent zones (i.e., Haroharo and Tarawera) and within the northern Ngakuru graben, and 2) we estimate an annual volume loss of 0.007 km3/yr, less than half the volume predicted by unconstrained inversions. Based on our best fitting model, we suggest that the observed seismicity in the vicinity of the caldera is driven by the long-term stressing from the contracting magma as well as periodic fluid and gas release into zones of increased permeability surrounding the cooling intrusion. Furthermore, we show that simple cooling models can explain long-term, near-linear, subsidence rates for periods of tens to hundreds of years.