  Extending recent results in ${\cal N}=2$ string compactifications, we propose
that the holomorphic anomaly equation satisfied by the modular completions of
the generating functions of refined BPS indices has a universal structure
independent of the number ${\cal N}$ of supersymmetries. We show that this
equation allows to recover all known results about modularity (under
$SL(2,\mathbb{Z})$ duality group) of BPS states in ${\cal N}=4$ string theory.
In particular, we reproduce the holomorphic anomaly characterizing the mock
modular behavior of quarter-BPS dyons and generalize it to the case of
non-trivial torsion invariant.
