clc
clear all;
format short g;
rho = 0.9;
t = 50; n = 20;
k = 100;
isim = 2000;

ic1 = 0; ic2 = 0;
hstr = randstream('mt19937ar', 'seed', 12345);  % set the seed number for random number generation

for is = 1:isim;
   u = randn(hstr,t+k,n);
   e = randn(hstr,t+k,n);
   
   y = u;
   x = e;
   
   for i = 2:t+k;
      y(i,:) = y(i-1,:)*rho + u(i,:);
      x(i,:) = x(i-1,:)*rho + e(i,:);
   end;
   
   y = y(k+1:t+k,:);
   x = x(k+1:t+k,:);
   
   y = y - repmat(mean(y),t,1);
   x = x - repmat(mean(x),t,1);
   
   xx = x(:);
   b = inv(xx'*xx)*xx'*y(:);
   uhat = y - x*b;
   ome1 = panelcov(x,uhat,4); %panel robust
   ome2 = panelcov(x,uhat,0); %ordinary t-ratio
    tr1 = b./sqrt(ome1);
    tr2 = b./sqrt(ome2);
    if abs(tr1) > 1.96; ic1 = ic1 + 1; end;
    if abs(tr2) > 1.96; ic2 = ic2 + 1; end;

end;
  [ ic1 ic2]
   
% break;

% Example of WRONG Bootstrap: starts ==========================
clear all;
  t = 50; k = 100; n = 20;
  rho = 0.9;
  % hstr = randstream('mt19937ar', 'seed', 12455);  % set the seed number for random number generation

   u = randn(t+k,n);
   e = randn(t+k,n);
   
   y = u;
   x = e;
   
   for i = 2:t+k;
      y(i,:) = y(i-1,:)*rho + u(i,:);
      x(i,:) = x(i-1,:)*rho + e(i,:);
   end;
   
   y = y(k+1:t+k,:);
   x = x(k+1:t+k,:);
   %y = y + x;
   dd = kron(eye(n),ones(t,1));
   
   xx = [dd x(:)];
   b = inv(xx'*xx)*xx'*y(:);
   db = repmat(b(1:n,:)',t,1);
   uhat = y - db - x*b(n+1);

   y1 = y - repmat(mean(y),t,1);
   x1 = x - repmat(mean(x),t,1);
   
   xx = x1(:);
   bb = inv(xx'*xx)*xx'*y1(:);
   uhat = y1 - x1*bb;
   ome1 = panelcov(x1,uhat,4); %panel robust
   ome2 = panelcov(x1,uhat,0); %ordinary t-ratio
    trat1 = bb./sqrt(ome1);
    trat2 = bb./sqrt(ome2);
    
   msim = 1000;
   csim = zeros(msim,2);
   
   for m = 1:msim;
       e = rand(t+k,1)*(t-1);
       e = 1 + floor(e);
       
       us = uhat(e,:);
       xs = x(e,:);
       us = us - repmat(mean(us),t+k,1); %recentering
       
       ys = us;
       for i = 2:t+k;
        %ys(i,:) = b(1:n,:)' + b(n+1)*xs(i,:) + us(i,:);
        ys(i,:) = b(1:n,:)' + us(i,:);
       end;
  
       ys = ys(k+1:t+k,:);
       xs = xs(k+1:t+k,:);
   
        ys = ys - repmat(mean(ys),t,1);
        xs = xs - repmat(mean(xs),t,1);
   
        xx = xs(:);
        bs = inv(xx'*xx)*xx'*ys(:);
        uh1 = ys - xs*bs;
        ome1 = panelcov(xs,uh1,4); %panel robust
        ome2 = panelcov(xs,uh1,0); %ordinary t-ratio
        tr1 = bs./sqrt(ome1);
        tr2 = bs./sqrt(ome2);
        csim(m,:) = [tr1 tr2];
   end;
        csim = sort(abs(csim));
        crtic = csim(0.95*msim,:);        
   
        [trat1 crtic(1)]
        [trat2 crtic(2)]