How to calculate time domain SNR using known sequence

I am using the following formula to calculate SNR of a real world complex baseband signal sampled at 1x Nyquist.

SNR      = Rxy(tm)^2 / [ Px*Py - Rxy(tm)^2 ]
SNR (dB) = 10*log10(SNR)

where

Rxy(tm) = peak of the cross correlation at time delay, tm 
Px = power in reference signal 
Py = power in received signal

I verified proper implementation of the formula using simulated real-valued and complex-valued signals with and without noise.

On real data the SNR estimates using the above formula are too low (by 10.0+ dB). I manually verified the actual SNR a few different ways. I used spectral analysis to visually measure the signal power to the noise floor. I also measured the signal power to noise power (when signal is off), and both of those techniques give me an answer closer to what I expect.

I am flummoxed as to why this equation is not working on real-world signals. Do I need to take the DC bias (mean of data) into account and add that back to the SNR estimate? If I do that then I get values closer to what I expect.

Reference: Formula came from Principles of Communications (Tranter, Ziemer) textbook

=Rxy(tm)^2 / (Py-Rxy(tm)^2)