A new window is opening in high-energy astronomy: X-ray polarimetry. With many missions currently under development and scheduled to launch as early as 2021, observations of the X-ray polarisation of accreting X-ray pulsars will soon be available. As polarisation is particularly sensitive to the geometry of the emission region, the upcoming polarimeters will shed new light on the emission mechanism of these objects, provided that we have sound theoretical models that agree with current spectroscopic and timing observation and that can make predictions of the polarisation parameters of the emission. We here present a new model for the polarised emission of accreting X-ray pulsars in the accretion column scenario that for the first time takes into account the macroscopic structure and dynamics of the accretion region and the propagation of the radiation toward the observer, including relativistic beaming, gravitational lensing and quantum electrodynamics. In this paper we present all the details of the model, while in a companion paper, we apply our model to predict the polarisation parameters of the bright X-ray pulsar Hercules X-1.