Extending recent results in mathcal{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 mathcal{N} of supersymmetries. We show that this equation allows to recover all known results about modularity (under SL(2, ℤ) duality group) of BPS states in mathcal{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.