Let $X$ be a complex $K3$ surface with an effective action of a group $G$ which preserves the holomorphic symplectic form. Let $$ Z_{X,G}(q) = \sum_{n=0}^{\infty} e\left(\operatorname{Hilb}^{n}(X)^{G} \right)\, q^{n-1} $$ be the generating function for the Euler characteristics of the Hilbert schemes of $G$-invariant length $n$ subschemes. We show that its reciprocal, $Z_{X,G}(q)^{-1}$ is the Fourier expansion of a modular cusp form of weight $\frac{1}{2} e(X/G)$ for the congruence subgroup $\Gamma_{0}(|G|)$. We give an explicit formula for $Z_{X,G}$ in terms of the Dedekind eta function for all 82 possible $(X,G)$. The key intermediate result we prove is of independent interest: it establishes an eta product identity for a certain shifted theta function of the root lattice of a simply laced root system. We extend our results to various refinements of the Euler characteristic, namely the Elliptic genus, the Chi-$y$ genus, and the motivic class.