function y = bbs2(myy, myx, num, k)
% non-overlapped block bootstrap
% num: times 
% k: block length
% assume a random walk for the DGP of myx
% 486 obs 

T = size(myy,1);
n = T/k; % the number of blocks
tboots = nan(num,1);

% % xcell1 = cell(T,1);
% % for ii = 1:T
% %     xcell1{ii} = eye(3);
% % end 
% % [behat,~,res,~] = mvregress(xcell1, myy, 'algorithm','mvn');

ymean = mean(myy,1);
res = myy - ymean;

% res1 = reshape(res(:,1),k,n);
% res2 = reshape(res(:,2),k,n);
% res3 = reshape(res(:,3),k,n);
res1 = reshape(res(:,1),k,n);
res2 = reshape(res(:,2),k,n);
res3 = reshape(res(:,3),k,n);
ypred = myy - res;
xinit1 = myx(1,1);
xinit2 = myx(1,2);
xinit3 = myx(1,3);
xdiff = [0 0 0; diff(myx,1)];
xdiff1 = reshape(xdiff(:,1),k,n);
xdiff2 = reshape(xdiff(:,2),k,n);
xdiff3 = reshape(xdiff(:,3),k,n);
xboots = nan(T,3);
xcell = cell(T,1);
% % ypred(:,1) = [ones(T,1),myx(:,1)]*[bhat(1); bhat(4)];
% % ypred(:,2) = [ones(T,1),myx(:,2)]*[bhat(2); bhat(4)];
% % ypred(:,3) = [ones(T,1),myx(:,3)]*[bhat(3); bhat(4)];
% 
% % xcell = cell(T,1);
% % for ii = 1:T
% %     xcell{ii} = [eye(3), myx(ii,:)'];
% % end

for i = 1:num
    ind1 = randi(n,1);
    ind = randi(n,n-1,1);
    xboots(:,1) = [xinit1;xdiff1(2:end,ind1);reshape(xdiff1(:,ind),[],1)];
    xboots(:,2) = [xinit2;xdiff2(2:end,ind1);reshape(xdiff2(:,ind),[],1)];
    xboots(:,3) = [xinit3;xdiff3(2:end,ind1);reshape(xdiff3(:,ind),[],1)];
    xboots = cumsum(xboots,1);
    for ii = 1:T
        xcell{ii} = [eye(3), xboots(ii,:)'];
    end
    myresid = [reshape(res1(:,[ind1;ind]),[],1),reshape(res2(:,[ind1;ind]),[],1),reshape(res3(:,[ind1;ind]),[],1)];
    yb = ypred + myresid;
%     xcell = cell(T,1);
%     for ii = 1:T
%         xcell{ii} = [eye(3), myx(ii,:)'];
%     end
    [b,~,~,CovB] = mvregress(xcell,yb, 'algorithm','mvn');
    tboots(i) = b(4)/sqrt(CovB(4,4));
    if mod(i,1000)==0
        disp(i)
    end
end

y = tboots;

end 