Power spectrum estimate from FFT
The difference is in the scaling of the power spectrum. I suppose you will find that the difference in scaling always equals $L/F_s$, which for the numbers in your question is indeed $8000/500=16$.

Take a look at Power Spectral Density Estimates Using FFT for the correct scaling. If you normalize the FFT result by the FFT length, you need to scale the squared magnitude of the normalized FFT by $L/F_s$. Furthermore, you need a factor of $2$ if you throw away the negative frequencies (I see you did that). Note, however, that you shouldn't apply a factor of $2$ for the bins at DC and at Nyquist, because they are not mirrored at negative frequencies.

Here is a simple Matlab code from the above quoted mathworks page for computing a periodogram-based one-sided power spectrum estimate using the FFT (my comments):


Fs = 1000;                          % sampling frequency (Hz)
N = length(x);                      % even! (=> bin N/2+1 is Nyquist)
xdft = fft(x);
xdft = xdft(1:N/2+1);               % DC to Nyquist (one-sided)
psdx = (1/(Fs*N)) * abs(xdft).^2;   % periodogram scaling
psdx(2:end-1) = 2*psdx(2:end-1);    % scaling for one-sided PSD
                                    % note: no factor 2 at DC and Nyquist