function part2simsamplingdistribution(K,L,N,sigma,numobs)

  %L is here assumed known
  
  c1 = N/K;
  c2 = L/N;
  c3 = L/K;
  
  angl = pi/4;
  actualmom2 = (2*pi/angl);
  actualmom3 = (2*pi/angl)^2;
  measuredi2 = [];
  measuredi3 = [];
  measuredalpha = [];
  obsvs = linspace(1,numobs,100);
  obsvs = floor(obsvs);
  for (obsnr=obsvs)
    obsnr
    samplecovsum1 = 0;
    samplecovsum2 = 0;
    samplecovsum3 = 0;
    for (thisobs=1:obsnr)
      V = zeros(N,L);
      for (k=1:L) 
        V(:,k) = ((exp(j*angl*rand(1))).^(0:(N-1)))';
      end
      V = (1/sqrt(N))*V;
      alpha = randn(L,K);
      nmatrix = sigma*randn(N,K);
      Y = V*alpha + nmatrix;
      W = (1/K)*Y*Y';%sample covariance matrix sample
      samplecovsum1 = samplecovsum1 + (1/N)*trace(W);
      samplecovsum2 = samplecovsum2 + (1/N)*trace(W^2);
      samplecovsum3 = samplecovsum3 + (1/N)*trace(W^3);
    end
    W1 = samplecovsum1/obsnr;
    W2 = samplecovsum2/obsnr;
    W3 = samplecovsum3/obsnr;
 
    I2 = ( W2     -      (c2*c3)...     
                  -      2*(sigma^2)*(c2+c3) ...    
                  -      (sigma^4)*(1+c1) ...
         )/(c2^2) - 1/c2;
    I3 = (W3      -      (3*(c2^2)*I2 + 3*c2*c3) ...
                  -      (3*(c2^2)*c3*I2 + c2*(c3^2) ) ...    
                  -      3*(sigma^2)*( (1+c1)*c2 ) ...    
                  -      3*(sigma^2)*( (1+c1)*(c2^2)*I2 + c3*(c3+2*c2) ) ...    
                  -      3*(sigma^4)*( c1^2+3*c1+1)*c2 ...     
                  -      (sigma^6)*(c1^2+3*c1+1) ...      
         )/(c2^3) - 1/(c2^2);
    measuredi2 = [measuredi2 I2];
    measuredi3 = [measuredi3 I3];
    measuredalpha = [measuredalpha  (2*pi/I2) ];
  end  
  
  plot(obsvs,measuredalpha,'k+',obsvs,angl*ones(1,length(obsvs)),'k-');
  xlabel('Number of observations');
  axis([0 numobs 0 1.4]);
  legend('Estimated \alpha','Actual \alpha');
  print -deps ../images/estimatedalpha.eps
  
  plot(obsvs,measuredi2,'k+',obsvs,actualmom2*ones(1,length(obsvs)),'k-',obsvs,measuredi3,'ko',obsvs,actualmom3*ones(1,length(obsvs)),'k--');
  xlabel('Number of observations');
  axis([0 numobs 0 120]);
  legend('Estimated I_2','Actual I_2','Estimated I_3','Actual I_3');
  print -deps ../images/estimatedi2i3.eps