Here are a couple of ways. Since the time function fills up the entire time window (as opposed to being some kind of pulse that dies off before reaching each end), it is a bit easier to get a close equality by using the rectangular rule approximation for each integral, [sum of points] x [width of interval]. The trapezoidal rule is not quite right for a periodic window.
You are using Fs for the frequency array step width, which conflicts a bit with the usual notation of Fs as the sampling frequency, but I kept Fs as you have it.
By introducing a time array and trapz, you are implicitly declaring the fft to be an approximate representation of a fourier transform integral. Accordingly in the expression for AA_S, there is a factor of Ts for the width of each interval in the riemann sum.
The TD and Freq results below agree.
N = 1000;
t = ((0:N-1)/N)*10;
Fs = 1/(N*Ts);
frequencyAxis = ((0:N-1)/N)*Fs;
Another way, which dispenses with interval widths and is more signal processing oriented, is
N = 1000;
r = rand(1,N);
E1 = sum(r.^2)
g = fft(r);
E2 = sum(abs(g).^2)/N
When you prove parseval's theorem and plug in ffts, there is a sum over the product of a couple of complex exponentials, and that sum is zero except for one instance where the product of the exponentials is 1. Then the sum over points gives N, which gets compensated for by the 1/N factor on the last llne. Needless to say if you reverse engineeer your way out of the first method above, you can obtain the second method.