clear;
clc;

t= 300; n=20;
k = 100;
u = randn(t+k,n);
 th = randn(t+k,1);
 u = repmat(th,1,n) + u;
y = u;
rho = 0.9;

for i = 2:t+k;
  y(i,:) = y(i-1,:)*rho + u(i,:);
end;

  y = y(k+1:t+k,:);
  
  y1 = y(2:t,:); y2 = y(1:t-1,:);
  y1 = y1 - repmat(mean(y1),t-1,1);
  y2 = y2 - repmat(mean(y2),t-1,1);
  
  b = sum(sum(y1.*y2))/sum(sum(y2.*y2));
  e = y1 - y2*b;
  e1 = e.*e; e1 = mean(mean(e1));
  tb = e1/sum(sum(y2.*y2));
  tb = (b-rho)./sqrt(tb);
  
  % =========== bootstrap ===================
  
  msim = 1000;
  csim = zeros(msim,1);
  
  for m = 1:msim;
       v = rand(t+k,1)*(t-1);
       v = 1 + floor(v);
       es = e(v,:);
   
       ys = es;
       for i = 2:t+k;
        ys(i,:) =  ys(i-1,:)*b + es(i,:);
       end;
  
       ys = ys(k+1:t+k,:);
   
       w1 = ys(2:t,:);
       w2 = ys(1:t-1,:);
   
       w1 = w1 - repmat(mean(w1),t-1,1);
       w2 = w2 - repmat(mean(w2),t-1,1);
  
        bw = sum(sum(w1.*w2))/sum(sum(w2.*w2));
        es = w1 - w2*bw;
        tbs = var(es(:))/sum(sum(w2.*w2));
        tbs = (bw-b)./sqrt(tbs);
  
       csim(m,:) = tbs;
   end;
  
   %bias = mean(csim) - b;
   %mue = b - bias;
   %[b mue bias]
   
   csim = sort(abs(csim));
   [tb csim(0.95*msim)]
   