We present an extension of the pair coupled cluster doubles (p-CCD) method to quasiparticles and apply it to the attractive pairing Hamiltonian. Near the transition point where number symmetry gets spontaneously broken, the proposed BCS-based p-CCD method yields significantly better energies than existing methods when compared to exact results obtained via solution of the Richardson equations. The quasiparticle p-CCD method has a low computational cost of O(N3) as a function of system size. This together with the high quality of results here demonstrated, points to considerable promise for the accurate description of strongly correlated systems with more realistic pairing interactions.