  Whenever available, refined BPS indices provide considerably more information
on the spectrum of BPS states than their unrefined version. Extending earlier
work on the modularity of generalized Donaldson-Thomas invariants counting
D4-D2-D0 brane bound states in type IIA strings on a Calabi-Yau threefold
$\mathfrak{Y}$, we construct the modular completion of generating functions of
refined BPS indices supported on a divisor class. Although for compact
$\mathfrak{Y}$ the refined indices are not protected, switching on the
refinement considerably simplifies the construction of the modular completion.
Furthermore, it leads to a non-commutative analogue of the TBA equations, which
suggests a quantization of the moduli space consistent with S-duality. In
contrast, for a local CY threefold given by the total space of the canonical
bundle over a complex surface $S$, refined BPS indices are well-defined, and
equal to Vafa-Witten invariants of $S$. Our construction provides a modular
completion of the generating function of these refined invariants for arbitrary
rank. In cases where all reducible components of the divisor class are
collinear (which occurs e.g. when $b_2(\mathfrak{Y})=1$, or in the local case),
we show that the holomorphic anomaly equation satisfied by the completed
generating function truncates at quadratic order. In the local case, it agrees
with an earlier proposal by Minahan et al for unrefined invariants, and extends
it to the refined level using the afore-mentioned non-commutative structure.
Finally, we show that these general predictions reproduce known results for
$U(2)$ and $U(3)$ Vafa-Witten theory on $\mathbb{P}^2$, and make them explicit
for $U(4)$.
