function part1vsim1(Llim,numavg)

Lvector = 5:5:Llim;

asympt = [ 1 2 5 44/3 ];

exactasympt = [];
momsexact = [];

for L=Lvector
  N = L
  exact(1) = 1;
  exact(2) = 2 -  1/N; 
  exact(3) = 5 - 6/N + 2/(N^2);
  exact(4) = 44/3 - 89/(3*N) + 133/(6*(N^2)) - 37/(6*(N^3));
  
  moms = zeros(1,4);
  for (l=1:numavg) 
    V = zeros(N,L);
    for (k=1:L) 
      V(:,k) = ((exp(2*pi*j*rand(1))).^(0:(N-1)))';
    end
    V = (1/sqrt(N))*V;
    moms(1) = moms(1) + (1/L)*trace((V'*V));
    moms(2) = moms(2) + (1/L)*trace((V'*V)^2);
    moms(3) = moms(3) + (1/L)*trace((V'*V)^3);
    moms(4) = moms(4) + (1/L)*trace((V'*V)^4);
  end
  moms = moms/numavg;
  
  exactasympt = [ exactasympt (exact - asympt) * (exact - asympt)' ];
  momsexact = [ momsexact (moms - exact) * (moms - exact)' ];
end

plot(Lvector,exactasympt,'k+',Lvector,momsexact,'ko');
legend_handle = legend('MSE between exact and asymptotic moments ','MSE between estimated and exact moments','Location','North');
set(legend_handle,'Box','off');
axis([5 Llim 0 1]);