In this paper, we study the Villain model in $\mathbb{Z}^d$ in dimension $d\geq 3$. It is conjectured, that the parallel correlation function in the infinite volume Gibbs state, i.e., the map $$ x \mapsto \langle \cosθ(0) \cosθ(x) \rangle_{μ_{\mathrm{Vil}, β}} -\left( \langle \cosθ(0) \rangle_{μ_{\mathrm{Vil}, β}} \right)^2, $$ decays like $|x|^{-2(d-2)}$ as $|x| \to \infty$ at low temperature. The results of Bricmont, Fontaine, Lebowitz, Lieb, and Spencer (1981) show that for the related XY model, this correlation decays at least as fast as $|x|^{2-d}$. We prove the optimal upper and lower bounds for the Villain model in $d=3$, up to a logarithmic correction, and also improve the upper bound in general dimensions. Our proof builds upon the approach developed in our previous article, which in turn is inspired by a key observation of Fröhlich and Spencer (1982): in the low temperature regime, a combination of duality transformation and renormalisation allows certain properties of the Villain model to be analysed in terms of a (vector-valued) $\nabla \varphi$ interface model. This latter model can be investigated using the Helffer-Sjöstrand representation formula combined with tools of elliptic and parabolic regularity.